Made compute snap fully parallel (11961084) · Commits · 郑智淋 / lammps

examples/snap/in.snap.compute

0 → 100644

+94 −0

Original line number	Diff line number	Diff line
		# Demonstrate bispectrum computes

		# Initialize simulation

		variable nsteps index 0
		variable nrep equal 1
		#variable a equal 3.316
		variable a equal 2.0
		units metal

		# generate the box and atom positions using a BCC lattice

		variable nx equal ${nrep}
		variable ny equal ${nrep}
		variable nz equal ${nrep}

		boundary p p p

		lattice bcc $a
		region box block 0 ${nx} 0 ${ny} 0 ${nz}
		create_box 2 box
		create_atoms 2 box

		mass * 180.88

		displace_atoms all random 0.1 0.1 0.1 123456

		# choose SNA parameters

		variable twojmax equal 2
		variable rcutfac equal 1.0
		variable rfac0 equal 0.99363
		variable rmin0 equal 0
		variable radelem1 equal 2.3
		variable radelem2 equal 2.0
		variable wj1 equal 1.0
		variable wj2 equal 0.96


		# Setup dummy potential to satisfy cutoff

		pair_style zero ${rcutfac}
		pair_coeff * *

		# Setup reference potential

		variable zblcutinner equal 4
		variable zblcutouter equal 4.8
		variable zblz equal 73
		pair_style zbl ${zblcutinner} ${zblcutouter}
		pair_coeff * * ${zblz} ${zblz}

		# set up old-style per-atom computes

		compute b all sna/atom ${rcutfac} ${rfac0} ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag 0 bzeroflag 0 switchflag 0
		compute vb all snav/atom ${rcutfac} ${rfac0} ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag 0 bzeroflag 0 switchflag 0
		compute db all snad/atom ${rcutfac} ${rfac0} ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag 0 bzeroflag 0 switchflag 0

		# perform sums over atoms

		group snapgroup1 type 1
		group snapgroup2 type 2
		compute bsum1 snapgroup1 reduce sum c_b[*]
		compute bsum2 snapgroup2 reduce sum c_b[*]
		# fix bsum1 all ave/time 1 1 1 c_bsum1 file bsum1.dat mode vector
		# fix bsum2 all ave/time 1 1 1 c_bsum2 file bsum2.dat mode vector
		compute vbsum all reduce sum c_vb[*]
		# fix vbsum all ave/time 1 1 1 c_vbsum file vbsum.dat mode vector

		# set up new-style global compute

		compute snap all snap ${rcutfac} ${rfac0} ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag 0 bzeroflag 0 switchflag 0
		compute snap_press all pressure NULL virial
		fix snap all ave/time 1 1 1 c_snap[*] file compute.snap.dat mode vector

		thermo 100

		# test output: 1: total potential energy
		# 2: xy component of stress tensor
		# 3: Sum(B_{000}^i, all i of type 2)
		# 4: xy component of Sum(Sum(r_j*dB_{222}^i/dr_j), all i of type 2), all j)
		# followed by counterparts from compute snap

		thermo_style custom &
		pe pxy c_bsum2[1] c_vbsum[60] &
		c_snap[1][11] c_snap[13][11] c_snap[1][6] c_snap[13][10]
		thermo_modify norm no

		# dump mydump_db all custom 1000 dump_db id c_db[*]
		# dump_modify mydump_db sort id

		# Run MD

		run ${nsteps}

src/SNAP/compute_snap.cpp

+49 −87

Original line number	Diff line number	Diff line
		@@ -45,7 +45,7 @@ enum{SCALAR,VECTOR,ARRAY};

		ComputeSnap::ComputeSnap(LAMMPS lmp, int narg, char *arg) :
		Compute(lmp, narg, arg), cutsq(NULL), list(NULL), snap(NULL),
		radelem(NULL), wjelem(NULL), snap_peratom(NULL)
		radelem(NULL), wjelem(NULL), snap_peratom(NULL), snapall(NULL)
		{

		array_flag = 1;
		@@ -136,9 +136,10 @@ ComputeSnap::ComputeSnap(LAMMPS lmp, int narg, char *arg) :
		natoms = atom->natoms;
		size_array_rows = 1+ndims_force*natoms+ndims_virial;
		size_array_cols = nperdim*atom->ntypes+1;
		lastcol = size_array_cols-1;

		ndims_peratom = ndims_force;
		size_peratom = ndims_peratomnperdimatom->ntypes;
		comm_reverse = size_peratom;

		nmax = 0;
		}
		@@ -148,6 +149,7 @@ ComputeSnap::ComputeSnap(LAMMPS lmp, int narg, char *arg) :
		ComputeSnap::~ComputeSnap()
		{
		memory->destroy(snap);
		memory->destroy(snapall);
		memory->destroy(snap_peratom);
		memory->destroy(radelem);
		memory->destroy(wjelem);
		@@ -185,7 +187,9 @@ void ComputeSnap::init()

		memory->create(snap,size_array_rows,size_array_cols,
		"snap:snap");
		array = snap;
		memory->create(snapall,size_array_rows,size_array_cols,
		"snap:snapall");
		array = snapall;

		// INCOMPLETE: modify->find_compute()
		// was called 223960 times by snappy Ta example
		@@ -235,10 +239,10 @@ void ComputeSnap::compute_array()
		}

		// clear global array
		// only need to zero out first row of bispectrum

		for (int icoeff = 0; icoeff < size_array_cols-1; icoeff++)
		snap[0][icoeff] = 0.0;
		for (int irow = 0; irow < size_array_rows; irow++)
		for (int icoeff = 0; icoeff < size_array_cols; icoeff++)
		snap[irow][icoeff] = 0.0;

		// clear local peratom array

		@@ -318,7 +322,7 @@ void ComputeSnap::compute_array()
		snaptr->rcutij[jj],jj);
		snaptr->compute_dbidrj();

		// Accumulate -dBi/dRi, -dBi/dRj
		// Accumulate dBi/dRi, -dBi/dRj

		double *snadi = snap_peratom[i]+typeoffset_local;
		double *snadj = snap_peratom[j]+typeoffset_local;
		@@ -404,91 +408,65 @@ void ComputeSnap::compute_array()
		}
		}

		// accumulate virial contributions before reverse compute

		dbdotr_compute();

		// communicate local peratom contributions between neighbor procs

		comm->reverse_comm_compute(this);

		// copy peratom data to global array
		// INCOMPLETE ignore local-global mapping for now
		// accumulate bispectrum force contributions to global array

		for (int itype = 0; itype < atom->ntypes; itype++) {
		const int typeoffset_local = ndims_peratomnperdimitype;
		const int typeoffset_global = nperdim*itype;
		for (int icoeff = 0; icoeff < nperdim; icoeff++) {
		int irow = 1;
		for (int i = 0; i < atom->nlocal; i++) {
		for (int i = 0; i < ntotal; i++) {
		double *snadi = snap_peratom[i]+typeoffset_local;
		snap[irow++][icoeff+typeoffset_global] = snadi[icoeff];
		snap[irow++][icoeff+typeoffset_global] = snadi[icoeff+yoffset];
		snap[irow++][icoeff+typeoffset_global] = snadi[icoeff+zoffset];
		int iglobal = atom->tag[i];
		int irow = 3*(iglobal-1)+1;
		snap[irow][icoeff+typeoffset_global] += snadi[icoeff];
		snap[irow+1][icoeff+typeoffset_global] += snadi[icoeff+yoffset];
		snap[irow+2][icoeff+typeoffset_global] += snadi[icoeff+zoffset];
		}
		}
		}

		// assign energy to last column

		int icol = size_array_cols-1;
		int irow = 0;
		double reference_energy = c_pe->compute_scalar();
		snap[irow++][icol] = reference_energy;

		// assign forces to last column
		// INCOMPLETE ignore local-global mapping for now
		// accumulate forces to global array

		for (int i = 0; i < atom->nlocal; i++) {
		snap[irow++][icol] = atom->f[i][0];
		snap[irow++][icol] = atom->f[i][1];
		snap[irow++][icol] = atom->f[i][2];
		int iglobal = atom->tag[i];
		int irow = 3*(iglobal-1)+1;
		snap[irow][lastcol] = atom->f[i][0];
		snap[irow+1][lastcol] = atom->f[i][1];
		snap[irow+2][lastcol] = atom->f[i][2];
		}

		// assign virial stress to last column
		// switch to Voigt notation
		// accumulate bispectrum virial contributions to global array

		c_virial->compute_vector();
		snap[irow++][icol] = c_virial->vector[0];
		snap[irow++][icol] = c_virial->vector[1];
		snap[irow++][icol] = c_virial->vector[2];
		snap[irow++][icol] = c_virial->vector[5];
		snap[irow++][icol] = c_virial->vector[4];
		snap[irow++][icol] = c_virial->vector[3];
		dbdotr_compute();

		}
		// sum up over all processes

		/* ---------------------------------------------------------------------- */
		MPI_Allreduce(&snap[0][0],&snapall[0][0],size_array_rows*size_array_cols,MPI_DOUBLE,MPI_SUM,world);

		int ComputeSnap::pack_reverse_comm(int n, int first, double *buf)
		{
		int i,m,last,icoeff;

		m = 0;
		last = first + n;
		for (i = first; i < last; i++)
		for (icoeff = 0; icoeff < size_peratom; icoeff++)
		buf[m++] = snap_peratom[i][icoeff];
		return m;
		}
		// assign energy to last column

		/* ---------------------------------------------------------------------- */
		int irow = 0;
		double reference_energy = c_pe->compute_scalar();
		snapall[irow++][lastcol] = reference_energy;

		void ComputeSnap::unpack_reverse_comm(int n, int list, double buf)
		{
		int i,j,m,icoeff;
		// assign virial stress to last column
		// switch to Voigt notation

		c_virial->compute_vector();
		irow += 3*natoms;
		snapall[irow++][lastcol] = c_virial->vector[0];
		snapall[irow++][lastcol] = c_virial->vector[1];
		snapall[irow++][lastcol] = c_virial->vector[2];
		snapall[irow++][lastcol] = c_virial->vector[5];
		snapall[irow++][lastcol] = c_virial->vector[4];
		snapall[irow++][lastcol] = c_virial->vector[3];

		m = 0;
		for (i = 0; i < n; i++) {
		j = list[i];
		for (icoeff = 0; icoeff < size_peratom; icoeff++)
		snap_peratom[j][icoeff] += buf[m++];
		}
		}

		/* ----------------------------------------------------------------------
		compute global virial contributions via summing dB dot r over
		own & ghost atoms, before reverse comm.
		compute global virial contributions via summing r_i.dB^j/dr_i over
		own & ghost atoms
		------------------------------------------------------------------------- */

		void ComputeSnap::dbdotr_compute()
		@@ -496,22 +474,8 @@ void ComputeSnap::dbdotr_compute()
		double **x = atom->x;
		int irow0 = 1+ndims_force*natoms;

		// zero virial entries

		for (int itype = 0; itype < atom->ntypes; itype++) {
		const int typeoffset_global = nperdim*itype;
		for (int icoeff = 0; icoeff < nperdim; icoeff++) {
		int irow = irow0;
		snap[irow++][icoeff+typeoffset_global] = 0.0;
		snap[irow++][icoeff+typeoffset_global] = 0.0;
		snap[irow++][icoeff+typeoffset_global] = 0.0;
		snap[irow++][icoeff+typeoffset_global] = 0.0;
		snap[irow++][icoeff+typeoffset_global] = 0.0;
		snap[irow++][icoeff+typeoffset_global] = 0.0;
		}
		}

		// sum over force on all particles including ghosts
		// sum over bispectrum contributions to forces
		// on all particles including ghosts

		int nall = atom->nlocal + atom->nghost;
		for (int i = 0; i < nall; i++)
		@@ -524,10 +488,6 @@ void ComputeSnap::dbdotr_compute()
		double dbdy = snadi[icoeff+yoffset];
		double dbdz = snadi[icoeff+zoffset];
		int irow = irow0;

		// RHS is xi*dSum(b_j, j in itype)/dxi
		// dSum(bj)/dxi is in first ntypes3nperdim elements of row snap_peratom[i]
		// LHS is Sum(RHS, i), each row of length ntypes*nperdim
		snap[irow++][icoeff+typeoffset_global] += dbdx*x[i][0];
		snap[irow++][icoeff+typeoffset_global] += dbdy*x[i][1];
		snap[irow++][icoeff+typeoffset_global] += dbdz*x[i][2];
		@@ -547,6 +507,8 @@ double ComputeSnap::memory_usage()

		double bytes = size_array_rowssize_array_cols
		sizeof(double); // snap
		bytes += size_array_rowssize_array_cols
		sizeof(double); // snapall
		bytes += nmaxsize_peratom sizeof(double); // snap_peratom
		bytes += snaptr->memory_usage(); // SNA object

src/SNAP/compute_snap.h

+2 −4

Original line number	Diff line number	Diff line
		@@ -31,17 +31,15 @@ class ComputeSnap : public Compute {
		void init();
		void init_list(int, class NeighList *);
		void compute_array();
		int pack_reverse_comm(int, int, double *);
		void unpack_reverse_comm(int, int , double );
		double memory_usage();

		private:
		int natoms, nmax, size_peratom;
		int natoms, nmax, size_peratom, lastcol;
		int ncoeff, nperdim, yoffset, zoffset;
		int ndims_peratom, ndims_force, ndims_virial;
		double **cutsq;
		class NeighList *list;
		double **snap;
		double snap, snapall;
		double **snap_peratom;
		double rcutfac;
		double *radelem;

Admin message