diff --git a/doc/Section_start.html b/doc/Section_start.html index cc9106ebf6..230c8644ed 100644 --- a/doc/Section_start.html +++ b/doc/Section_start.html @@ -189,6 +189,7 @@ within the LAMMPS code. The options that are currently recogized are:

The read_data and dump commands will read/write gzipped files if you -compile with -DLAMMPS_GZIP. It requires that your Unix support the -"popen" command. +compile with -DLAMMPS_GZIP. It requires that your machine supports +the "popen" function in the standard runtime library and that a gzip +executable can be found by LAMMPS during a run.

-

If you use -DLAMMPS_JPEG, the dump image command will be -able to write out JPEG image files. For JPEG files, you must also link -LAMMPS with a JPEG library, as described below. If you use +

If you use -DLAMMPS_JPEG, the dump image command +will be able to write out JPEG image files. For JPEG files, you must +also link LAMMPS with a JPEG library, as described below. If you use -DLAMMPS_PNG, the dump image command will be able to write out PNG image files. For PNG files, you must also link LAMMPS with a PNG library, as described below. If neither of those two defines are -used, LAMMPS will only be able to write out text-based PPM image +used, LAMMPS will only be able to write out uncompressed PPM image files.

+

If you use -DLAMMPS_FFMPEG, the dump movie command +will be available to support on-the-fly generation of rendered movies +the need to store intermediate image files. It requires that your +machines supports the "popen" function in the standard runtime library +and that an FFmpeg executable can be found by LAMMPS during the run. +

Using -DLAMMPS_MEMALIGN= enables the use of the posix_memalign() call instead of malloc() when large chunks or memory are allocated by LAMMPS. This can help to make more efficient use of diff --git a/doc/Section_start.txt b/doc/Section_start.txt index f905bd2e2a..f785d4d103 100644 --- a/doc/Section_start.txt +++ b/doc/Section_start.txt @@ -183,6 +183,7 @@ within the LAMMPS code. The options that are currently recogized are: -DLAMMPS_GZIP -DLAMMPS_JPEG -DLAMMPS_PNG +-DLAMMPS_FFMPEG -DLAMMPS_MEMALIGN -DLAMMPS_XDR -DLAMMPS_SMALLBIG @@ -194,18 +195,25 @@ within the LAMMPS code. The options that are currently recogized are: -DPACK_MEMCPY :ul The read_data and dump commands will read/write gzipped files if you -compile with -DLAMMPS_GZIP. It requires that your Unix support the -"popen" command. +compile with -DLAMMPS_GZIP. It requires that your machine supports +the "popen" function in the standard runtime library and that a gzip +executable can be found by LAMMPS during a run. -If you use -DLAMMPS_JPEG, the "dump image"_dump.html command will be -able to write out JPEG image files. For JPEG files, you must also link -LAMMPS with a JPEG library, as described below. If you use +If you use -DLAMMPS_JPEG, the "dump image"_dump_image.html command +will be able to write out JPEG image files. For JPEG files, you must +also link LAMMPS with a JPEG library, as described below. If you use -DLAMMPS_PNG, the "dump image"_dump.html command will be able to write out PNG image files. For PNG files, you must also link LAMMPS with a PNG library, as described below. If neither of those two defines are -used, LAMMPS will only be able to write out text-based PPM image +used, LAMMPS will only be able to write out uncompressed PPM image files. +If you use -DLAMMPS_FFMPEG, the "dump movie"_dump_image.html command +will be available to support on-the-fly generation of rendered movies +the need to store intermediate image files. It requires that your +machines supports the "popen" function in the standard runtime library +and that an FFmpeg executable can be found by LAMMPS during the run. + Using -DLAMMPS_MEMALIGN= enables the use of the posix_memalign() call instead of malloc() when large chunks or memory are allocated by LAMMPS. This can help to make more efficient use of diff --git a/doc/dump.html b/doc/dump.html index eb83e900ad..3d2e4f4a24 100644 --- a/doc/dump.html +++ b/doc/dump.html @@ -97,8 +97,8 @@ dump 2 subgroup atom 50 dump.run.bin dump 4a all custom 100 dump.myforce.* id type x y vx fx dump 4b flow custom 100 dump.%.myforce id type c_myF[3] v_ke -dump 2 inner cfg 10 dump.snap.*.cfg id type xs ys zs vx vy vz -dump snap all cfg 100 dump.config.*.cfg id type xs ys zs id type c_Stress2 +dump 2 inner cfg 10 dump.snap.*.cfg mass type xs ys zs vx vy vz +dump snap all cfg 100 dump.config.*.cfg mass type xs ys zs id type c_Stress2 dump 1 all xtc 1000 file.xtc dump e_data all custom 100 dump.eff id type x y z spin eradius fx fy fz eforce @@ -237,9 +237,9 @@ extended CFG format files, as used by the package. Since the extended CFG format uses a single snapshot of the system per file, a wildcard "*" must be included in the filename, as discussed below. The list of atom attributes for style cfg must -begin with either "id type xs ys zs" or "id type xsu ysu zsu" or +begin with either "mass type xs ys zs" or "mass type xsu ysu zsu" since these quantities are needed to -write the CFG files in the appropriate format (though the "id" and +write the CFG files in the appropriate format (though the "mass" and "type" fields do not appear explicitly in the file). Any remaining attributes will be stored as "auxiliary properties" in the CFG files. Note that you will typically want to use the dump_modify diff --git a/doc/dump.txt b/doc/dump.txt index 743acf148a..f68415f422 100644 --- a/doc/dump.txt +++ b/doc/dump.txt @@ -85,8 +85,8 @@ dump myDump all atom 100 dump.atom dump 2 subgroup atom 50 dump.run.bin dump 4a all custom 100 dump.myforce.* id type x y vx fx dump 4b flow custom 100 dump.%.myforce id type c_myF\[3\] v_ke -dump 2 inner cfg 10 dump.snap.*.cfg id type xs ys zs vx vy vz -dump snap all cfg 100 dump.config.*.cfg id type xs ys zs id type c_Stress[2] +dump 2 inner cfg 10 dump.snap.*.cfg mass type xs ys zs vx vy vz +dump snap all cfg 100 dump.config.*.cfg mass type xs ys zs id type c_Stress[2] dump 1 all xtc 1000 file.xtc dump e_data all custom 100 dump.eff id type x y z spin eradius fx fy fz eforce :pre @@ -225,9 +225,9 @@ extended CFG format files, as used by the package. Since the extended CFG format uses a single snapshot of the system per file, a wildcard "*" must be included in the filename, as discussed below. The list of atom attributes for style {cfg} must -begin with either "id type xs ys zs" or "id type xsu ysu zsu" or +begin with either "mass type xs ys zs" or "mass type xsu ysu zsu" since these quantities are needed to -write the CFG files in the appropriate format (though the "id" and +write the CFG files in the appropriate format (though the "mass" and "type" fields do not appear explicitly in the file). Any remaining attributes will be stored as "auxiliary properties" in the CFG files. Note that you will typically want to use the "dump_modify diff --git a/doc/fix_phonon.html b/doc/fix_phonon.html index 0c46909990..e013e09928 100644 --- a/doc/fix_phonon.html +++ b/doc/fix_phonon.html @@ -25,8 +25,14 @@

  • Nwait = wait this many timesteps before measuring -
  • map_file = file containing the mapping info between atom ID and the lattice indices +
  • map_file = file or GAMMA +
      file is the file that contains the mapping info between atom ID and the lattice indices. 
+
    +
      GAMMA flags to treate the whole simulation box as a unit cell, so that the mapping
+  info can be generated internally. In this case, dynamical matrix at only the gamma-point
+  will/can be evaluated. 
+
  • prefix = prefix for output files
  • one or none keyword/value pairs may be appended @@ -43,7 +49,8 @@

    Examples: 

    fix 1 all phonon 20 5000 200000 map.in LJ1D sysdim 1
-fix 1 all phonon 20 5000 200000 map.in EAM3D 
+fix 1 all phonon 20 5000 200000 map.in EAM3D
+fix 1 all phonon 10 5000 500000 GAMMA EAM0D nasr 100

    Description:

    @@ -124,6 +131,13 @@ which lattice point; the lattice indices start from 0. An auxiliary code,     latgen, can be employed to generate the compatible map file for various crystals.

    +

    In case one simulates an aperiodic system, where the whole simulation box +is treated as a unit cell, one can set map_file as GAMMA, so that the mapping +info will be generated internally and a file is not needed. In this case, the +dynamical matrix at only the gamma-point will/can be evaluated. Please keep in +mind that fix-phonon is designed for cyrstals, it will be inefficient and +even degrade the performance of lammps in case the unit cell is too large. +

    The calculated dynamical matrix elements are written out in energy/distance^2/mass units. The coordinates for q points in the log file is in the units of the basis vectors of the diff --git a/doc/fix_phonon.txt b/doc/fix_phonon.txt index 9590ded071..50285d06be 100644 --- a/doc/fix_phonon.txt +++ b/doc/fix_phonon.txt @@ -17,7 +17,12 @@ phonon = style name of this fix command :l N = measure the Green's function every this many timesteps :l Noutput = output the dynamical matrix every this many measurements :l Nwait = wait this many timesteps before measuring :l -map_file = file containing the mapping info between atom ID and the lattice indices :l +map_file = {file} or {GAMMA} :l + {file} is the file that contains the mapping info between atom ID and the lattice indices. :pre + + {GAMMA} flags to treate the whole simulation box as a unit cell, so that the mapping + info can be generated internally. In this case, dynamical matrix at only the gamma-point + will/can be evaluated. :pre prefix = prefix for output files :l one or none keyword/value pairs may be appended :l keyword = {sysdim} or {nasr} :l @@ -30,7 +35,8 @@ keyword = {sysdim} or {nasr} :l [Examples:] fix 1 all phonon 20 5000 200000 map.in LJ1D sysdim 1 -fix 1 all phonon 20 5000 200000 map.in EAM3D :pre +fix 1 all phonon 20 5000 200000 map.in EAM3D +fix 1 all phonon 10 5000 500000 GAMMA EAM0D nasr 100 :pre [Description:] @@ -111,6 +117,13 @@ which lattice point; the lattice indices start from 0. An auxiliary code, "latgen"_http://code.google.com/p/latgen, can be employed to generate the compatible map file for various crystals. +In case one simulates an aperiodic system, where the whole simulation box +is treated as a unit cell, one can set {map_file} as {GAMMA}, so that the mapping +info will be generated internally and a file is not needed. In this case, the +dynamical matrix at only the gamma-point will/can be evaluated. Please keep in +mind that fix-phonon is designed for cyrstals, it will be inefficient and +even degrade the performance of lammps in case the unit cell is too large. + The calculated dynamical matrix elements are written out in "energy/distance^2/mass"_units.html units. The coordinates for {q} points in the log file is in the units of the basis vectors of the diff --git a/src/USER-PHONON/fix_phonon.cpp b/src/USER-PHONON/fix_phonon.cpp index 83f6679076..5e3113300a 100644 --- a/src/USER-PHONON/fix_phonon.cpp +++ b/src/USER-PHONON/fix_phonon.cpp @@ -1,927 +1,927 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ -/* ---------------------------------------------------------------------- - Contributing authors: - Ling-Ti Kong - - Contact: - School of Materials Science and Engineering, - Shanghai Jiao Tong University, - 800 Dongchuan Road, Minhang, - Shanghai 200240, CHINA - - konglt@sjtu.edu.cn; konglt@gmail.com -------------------------------------------------------------------------- */ - -#include "string.h" -#include "fix_phonon.h" -#include "atom.h" -#include "compute.h" -#include "domain.h" -#include "fft3d_wrap.h" -#include "force.h" -#include "group.h" -#include "lattice.h" -#include "modify.h" -#include "update.h" -#include "citeme.h" -#include "memory.h" -#include "error.h" - -using namespace LAMMPS_NS; -using namespace FixConst; - -#define INVOKED_SCALAR 1 -#define INVOKED_VECTOR 2 -#define MAXLINE 512 - -static const char cite_fix_phonon[] = - "fix phonon command:\n\n" - "@Article{Kong11,\n" - " author = {L. T. Kong},\n" - " title = {Phonon dispersion measured directly from molecular dynamics simulations},\n" - " journal = {Comp.~Phys.~Comm.},\n" - " year = 2011,\n" - " volume = 182,\n" - " pages = {2201--2207}\n" - "}\n\n"; - -/* ---------------------------------------------------------------------- */ - -FixPhonon::FixPhonon(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) -{ - if (lmp->citeme) lmp->citeme->add(cite_fix_phonon); - - MPI_Comm_rank(world,&me); - MPI_Comm_size(world,&nprocs); - - if (narg < 8) error->all(FLERR,"Illegal fix phonon command: number of arguments < 8"); - - nevery = force->inumeric(FLERR, arg[3]); // Calculate this fix every n steps! - if (nevery < 1) error->all(FLERR,"Illegal fix phonon command"); - - nfreq = force->inumeric(FLERR, arg[4]); // frequency to output result - if (nfreq < 1) error->all(FLERR,"Illegal fix phonon command"); - - waitsteps = ATOBIGINT(arg[5]); // Wait this many timesteps before actually measuring - if (waitsteps < 0) error->all(FLERR,"Illegal fix phonon command: waitsteps < 0 !"); - - int n = strlen(arg[6]) + 1; // map file - mapfile = new char[n]; - strcpy(mapfile, arg[6]); - - n = strlen(arg[7]) + 1; // prefix of output - prefix = new char[n]; - strcpy(prefix, arg[7]); - logfile = new char[n+4]; - sprintf(logfile,"%s.log",prefix); - - int sdim = sysdim = domain->dimension; - int iarg = 8; - nasr = 20; - - // other command line options - while (iarg < narg){ - if (strcmp(arg[iarg],"sysdim") == 0){ - if (++iarg >= narg) error->all(FLERR,"Illegal fix phonon command: incomplete command line options."); - sdim = force->inumeric(FLERR, arg[iarg]); - if (sdim < 1) error->all(FLERR,"Illegal fix phonon command: sysdim should not be less than 1."); - - } else if (strcmp(arg[iarg],"nasr") == 0){ - if (++iarg >= narg) error->all(FLERR,"Illegal fix phonon command: incomplete command line options."); - nasr = force->inumeric(FLERR, arg[iarg]); - - } else { - error->all(FLERR,"Illegal fix phonon command: unknown option read!"); - } - - ++iarg; - } - - // get the dimension of the simulation; 1D is possible by specifying the option of "sysdim 1" - if (sdim < sysdim) sysdim = sdim; - nasr = MAX(0, nasr); - - // get the total number of atoms in group - ngroup = static_cast(group->count(igroup)); - if (ngroup < 1) error->all(FLERR,"No atom found for fix phonon!"); - - // MPI gatherv related variables - recvcnts = new int[nprocs]; - displs = new int[nprocs]; - - // mapping index - tag2surf.clear(); // clear map info - surf2tag.clear(); - - // get the mapping between lattice indices and atom IDs - readmap(); delete []mapfile; - if (nucell == 1) nasr = MIN(1,nasr); - - // get the mass matrix for dynamic matrix - getmass(); - - // create FFT and allocate memory for FFT - // here the parallization is done on the x direction only - nxlo = 0; - int *nx_loc = new int [nprocs]; - for (int i = 0; i < nprocs; ++i){ - nx_loc[i] = nx/nprocs; - if (i < nx%nprocs) ++nx_loc[i]; - } - for (int i = 0; i < me; ++i) nxlo += nx_loc[i]; - nxhi = nxlo + nx_loc[me] - 1; - mynpt = nx_loc[me]*ny*nz; - mynq = mynpt; - - fft_dim = nucell*sysdim; - fft_dim2 = fft_dim*fft_dim; - fft_nsend = mynpt*fft_dim; - - fft_cnts = new int[nprocs]; - fft_disp = new int[nprocs]; - fft_disp[0] = 0; - for (int i = 0; i < nprocs; ++i) fft_cnts[i] = nx_loc[i]*ny*nz*fft_dim; - for (int i = 1; i < nprocs; ++i) fft_disp[i] = fft_disp[i-1]+fft_cnts[i-1]; - delete []nx_loc; - - fft = new FFT3d(lmp,world,nz,ny,nx,0,nz-1,0,ny-1,nxlo,nxhi,0,nz-1,0,ny-1,nxlo,nxhi,0,0,&mysize); - memory->create(fft_data, MAX(1,mynq)*2, "fix_phonon:fft_data"); - - // allocate variables; MAX(1,... is used because NULL buffer will result in error for MPI - memory->create(RIloc,ngroup,(sysdim+1),"fix_phonon:RIloc"); - memory->create(RIall,ngroup,(sysdim+1),"fix_phonon:RIall"); - memory->create(Rsort,ngroup, sysdim, "fix_phonon:Rsort"); - - memory->create(Rnow, MAX(1,mynpt),fft_dim,"fix_phonon:Rnow"); - memory->create(Rsum, MAX(1,mynpt),fft_dim,"fix_phonon:Rsum"); - - memory->create(basis,nucell, sysdim, "fix_phonon:basis"); - - // because of hermit, only nearly half of q points are stored - memory->create(Rqnow,MAX(1,mynq),fft_dim, "fix_phonon:Rqnow"); - memory->create(Rqsum,MAX(1,mynq),fft_dim2,"fix_phonon:Rqsum"); - memory->create(Phi_q,MAX(1,mynq),fft_dim2,"fix_phonon:Phi_q"); - - // variable to collect all local Phi to root - if (me == 0) memory->create(Phi_all,ntotal,fft_dim2,"fix_phonon:Phi_all"); - else memory->create(Phi_all,1,1,"fix_phonon:Phi_all"); - - // output some information on the system to log file - if (me == 0){ - flog = fopen(logfile, "w"); - if (flog == NULL) { - char str[MAXLINE]; - sprintf(str,"Can not open output file %s",logfile); - error->one(FLERR,str); - } - for (int i = 0; i < 60; ++i) fprintf(flog,"#"); fprintf(flog,"\n"); - fprintf(flog,"# group name of the atoms under study : %s\n", group->names[igroup]); - fprintf(flog,"# total number of atoms in the group : %d\n", ngroup); - fprintf(flog,"# dimension of the system : %d D\n", sysdim); - fprintf(flog,"# number of atoms per unit cell : %d\n", nucell); - fprintf(flog,"# dimension of the FFT mesh : %d x %d x %d\n", nx, ny, nz); - fprintf(flog,"# number of wait steps before measurement : " BIGINT_FORMAT "\n", waitsteps); - fprintf(flog,"# frequency of the measurement : %d\n", nevery); - fprintf(flog,"# output result after this many measurement: %d\n", nfreq); - fprintf(flog,"# number of processors used by this run : %d\n", nprocs); - for (int i = 0; i < 60; ++i) fprintf(flog,"#"); fprintf(flog,"\n"); - fprintf(flog,"# mapping information between lattice index and atom id\n"); - fprintf(flog,"# nx ny nz nucell\n"); - fprintf(flog,"%d %d %d %d\n", nx, ny, nz, nucell); - fprintf(flog,"# l1 l2 l3 k atom_id\n"); - int ix, iy, iz, iu; - for (idx = 0; idx < ngroup; ++idx){ - itag = surf2tag[idx]; - iu = idx%nucell; - iz = (idx/nucell)%nz; - iy = (idx/(nucell*nz))%ny; - ix = (idx/(nucell*nz*ny))%nx; - fprintf(flog,"%d %d %d %d %d\n", ix, iy, iz, iu, itag); - } - for (int i = 0; i < 60; ++i) fprintf(flog,"#"); fprintf(flog,"\n"); - fflush(flog); - } - surf2tag.clear(); - - // default temperature is from thermo - TempSum = new double[sysdim]; - id_temp = new char[12]; - strcpy(id_temp,"thermo_temp"); - int icompute = modify->find_compute(id_temp); - temperature = modify->compute[icompute]; - inv_nTemp = 1.0/group->count(temperature->igroup); - -return; -} // end of constructor - -/* ---------------------------------------------------------------------- */ - -void FixPhonon::post_run() -{ - // compute and output final results - if (ifreq > 0 && ifreq != nfreq) postprocess(); - if (me == 0) fclose(flog); -} - -/* ---------------------------------------------------------------------- */ - -FixPhonon::~FixPhonon() -{ - // delete locally stored array - memory->destroy(RIloc); - memory->destroy(RIall); - memory->destroy(Rsort); - memory->destroy(Rnow); - memory->destroy(Rsum); - - memory->destroy(basis); - - memory->destroy(Rqnow); - memory->destroy(Rqsum); - memory->destroy(Phi_q); - memory->destroy(Phi_all); - - delete []recvcnts; - delete []displs; - delete []prefix; - delete []logfile; - delete []fft_cnts; - delete []fft_disp; - delete []id_temp; - delete []TempSum; - delete []M_inv_sqrt; - delete []basetype; - - // destroy FFT - delete fft; - memory->sfree(fft_data); - - // clear map info - tag2surf.clear(); - surf2tag.clear(); - -} - -/* ---------------------------------------------------------------------- */ - -int FixPhonon::setmask() -{ - int mask = 0; - mask |= END_OF_STEP; - -return mask; -} - -/* ---------------------------------------------------------------------- */ - -void FixPhonon::init() -{ - // warn if more than one fix-phonon - int count = 0; - for (int i = 0; i < modify->nfix; ++i) if (strcmp(modify->fix[i]->style,"gfc") == 0) ++count; - if (count > 1 && me == 0) error->warning(FLERR,"More than one fix phonon defined"); // just warn, but allowed. - -return; -} - -/* ---------------------------------------------------------------------- */ - -void FixPhonon::setup(int flag) -{ - // initialize accumulating variables - for (int i = 0; i < sysdim; ++i) TempSum[i] = 0.; - - for (int i = 0; i < mynpt; ++i) - for (int j = 0; j < fft_dim; ++j) Rsum[i][j] = 0.; - - for (int i =0; i < mynq; ++i) - for (int j =0; j < fft_dim2; ++j) Rqsum[i][j] = std::complex (0.,0.); - - for (int i = 0; i < 6; ++i) hsum[i] = 0.; - - for (int i = 0; i < nucell; ++i) - for (int j = 0; j < sysdim; ++j) basis[i][j] = 0.; - - prev_nstep = update->ntimestep; - neval = ifreq = 0; - -return; -} - -/* ---------------------------------------------------------------------- */ - -void FixPhonon::end_of_step() -{ - if ( (update->ntimestep-prev_nstep) <= waitsteps) return; - - double **x = atom->x; - int *mask = atom->mask; - int *tag = atom->tag; - int *image = atom->image; - int nlocal = atom->nlocal; - - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; - double *h = domain->h; - double xbox, ybox, zbox; - - int i,idim,jdim,ndim; - double xcur[3]; - - // to get the current temperature - if (!(temperature->invoked_flag & INVOKED_VECTOR)) temperature->compute_vector(); - for (idim = 0; idim < sysdim; ++idim) TempSum[idim] += temperature->vector[idim]; - - // evaluate R(r) on local proc - nfind = 0; - for (i = 0; i < nlocal; ++i){ - if (mask[i] & groupbit){ - itag = tag[i]; - idx = tag2surf[itag]; - - domain->unmap(x[i], image[i], xcur); - - for (idim = 0; idim < sysdim; ++idim) RIloc[nfind][idim] = xcur[idim]; - RIloc[nfind++][sysdim] = double(idx); - } - } - - // gather R(r) on local proc, then sort and redistribute to all procs for FFT - nfind *= (sysdim+1); - displs[0] = 0; - for (i = 0; i < nprocs; ++i) recvcnts[i] = 0; - MPI_Gather(&nfind,1,MPI_INT,recvcnts,1,MPI_INT,0,world); - for (i = 1; i < nprocs; ++i) displs[i] = displs[i-1] + recvcnts[i-1]; - - MPI_Gatherv(RIloc[0],nfind,MPI_DOUBLE,RIall[0],recvcnts,displs,MPI_DOUBLE,0,world); - if (me == 0){ - for (i = 0; i < ngroup; ++i){ - idx = static_cast(RIall[i][sysdim]); - for (idim = 0; idim < sysdim; ++idim) Rsort[idx][idim] = RIall[i][idim]; - } - } - MPI_Scatterv(Rsort[0],fft_cnts,fft_disp, MPI_DOUBLE, Rnow[0], fft_nsend, MPI_DOUBLE,0,world); - - // get Rsum - for (idx = 0; idx < mynpt; ++idx) - for (idim = 0; idim < fft_dim; ++idim) Rsum[idx][idim] += Rnow[idx][idim]; - - // FFT R(r) to get R(q) - for (idim = 0; idim < fft_dim; ++idim){ - int m=0; - for (idx = 0; idx < mynpt; ++idx){ - fft_data[m++] = Rnow[idx][idim]; - fft_data[m++] = 0.; - } - - fft->compute(fft_data, fft_data, -1); - - m = 0; - for (idq = 0; idq < mynq; ++idq){ - Rqnow[idq][idim] = std::complex(fft_data[m], fft_data[m+1]); - m += 2; - } - } - - // to get sum(R(q).R(q)*) - for (idq = 0; idq < mynq; ++idq){ - ndim = 0; - for (idim = 0; idim < fft_dim; ++idim) - for (jdim = 0; jdim < fft_dim; ++jdim) Rqsum[idq][ndim++] += Rqnow[idq][idim]*conj(Rqnow[idq][jdim]); - } - - // get basis info - if (fft_dim > sysdim){ - double dist2orig[3]; - for (idx = 0; idx < mynpt; ++idx){ - ndim = sysdim; - for (i = 1; i < nucell; ++i){ - for (idim = 0; idim < sysdim; ++idim) dist2orig[idim] = Rnow[idx][ndim++] - Rnow[idx][idim]; - domain->minimum_image(dist2orig); - for (idim = 0; idim < sysdim; ++idim) basis[i][idim] += dist2orig[idim]; - } - } - } - // get lattice vector info - for (int i = 0; i < 6; ++i) hsum[i] += h[i]; - - // increment counter - ++neval; - - // compute and output Phi_q after every nfreq evaluations - if (++ifreq == nfreq) postprocess(); - -return; -} // end of end_of_step() - -/* ---------------------------------------------------------------------- */ - -double FixPhonon::memory_usage() -{ - double bytes = sizeof(double)*2*mynq - + sizeof(std::map)*2*ngroup - + sizeof(double)*(ngroup*(3*sysdim+2)+mynpt*fft_dim*2) - + sizeof(std::complex)*MAX(1,mynq)*fft_dim *(1+2*fft_dim) - + sizeof(std::complex)*ntotal*fft_dim2 - + sizeof(int) * nprocs * 4; - return bytes; -} - -/* ---------------------------------------------------------------------- */ - -int FixPhonon::modify_param(int narg, char **arg) -{ - if (strcmp(arg[0],"temp") == 0) { - if (narg < 2) error->all(FLERR,"Illegal fix_modify command"); - delete [] id_temp; - int n = strlen(arg[1]) + 1; - id_temp = new char[n]; - strcpy(id_temp,arg[1]); - - int icompute = modify->find_compute(id_temp); - if (icompute < 0) error->all(FLERR,"Could not find fix_modify temp ID"); - temperature = modify->compute[icompute]; - - if (temperature->tempflag == 0) - error->all(FLERR,"Fix_modify temp ID does not compute temperature"); - inv_nTemp = 1.0/group->count(temperature->igroup); - - return 2; - } - -return 0; -} - -/* ---------------------------------------------------------------------- - * private method, to get the mass matrix for dynamic matrix - * --------------------------------------------------------------------*/ -void FixPhonon::getmass() -{ - int nlocal = atom->nlocal; - int *mask = atom->mask; - int *tag = atom->tag; - int *type = atom->type; - double *rmass = atom->rmass; - double *mass = atom->mass; - double *mass_one, *mass_all; - double *type_one, *type_all; - - mass_one = new double[nucell]; - mass_all = new double[nucell]; - type_one = new double[nucell]; - type_all = new double[nucell]; - for (int i = 0; i < nucell; ++i) mass_one[i] = type_one[i] = 0.; - - if (rmass){ - for (int i = 0; i < nlocal; ++i){ - if (mask[i] & groupbit){ - itag = tag[i]; - idx = tag2surf[itag]; - int iu = idx%nucell; - mass_one[iu] += rmass[i]; - type_one[iu] += double(type[i]); - } - } - } else { - for (int i = 0; i < nlocal; ++i){ - if (mask[i] & groupbit){ - itag = tag[i]; - idx = tag2surf[itag]; - int iu = idx%nucell; - mass_one[iu] += mass[type[i]]; - type_one[iu] += double(type[i]); - } - } - } - - MPI_Allreduce(mass_one,mass_all,nucell,MPI_DOUBLE,MPI_SUM,world); - MPI_Allreduce(type_one,type_all,nucell,MPI_DOUBLE,MPI_SUM,world); - - M_inv_sqrt = new double[nucell]; - basetype = new int[nucell]; - - double inv_total = 1./double(ntotal); - for (int i = 0; i < nucell; ++i){ - mass_all[i] *= inv_total; - M_inv_sqrt[i] = sqrt(1./mass_all[i]); - - basetype[i] = int(type_all[i]*inv_total); - } - delete []mass_one; - delete []mass_all; - delete []type_one; - delete []type_all; - -return; -} - - -/* ---------------------------------------------------------------------- - * private method, to read the mapping info from file - * --------------------------------------------------------------------*/ - -void FixPhonon::readmap() -{ - int info = 0; - - // auto-generate mapfile for "cluster" (gamma only system) - if (strcmp(mapfile, "GAMMA") == 0){ - nx = ny = nz = ntotal = 1; - nucell = ngroup; - - int *tag_loc, *tag_all; - memory->create(tag_loc,ngroup,"fix_phonon:tag_loc"); - memory->create(tag_all,ngroup,"fix_phonon:tag_all"); - - // get atom IDs on local proc - int nfind = 0; - for (int i = 0; i < atom->nlocal; ++i){ - if (atom->mask[i] & groupbit) tag_loc[nfind++] = atom->tag[i]; - } - - // gather IDs on local proc - displs[0] = 0; - for (int i = 0; i < nprocs; ++i) recvcnts[i] = 0; - MPI_Allgather(&nfind,1,MPI_INT,recvcnts,1,MPI_INT,world); - for (int i = 1; i < nprocs; ++i) displs[i] = displs[i-1] + recvcnts[i-1]; - - MPI_Allgatherv(tag_loc,nfind,MPI_INT,tag_all,recvcnts,displs,MPI_INT,world); - for (int i = 0; i < ngroup; ++i){ - itag = tag_all[i]; - tag2surf[itag] = i; - surf2tag[i] = itag; - } - - memory->destroy(tag_loc); - memory->destroy(tag_all); - return; - } - - // read from map file for others - char line[MAXLINE]; - FILE *fp = fopen(mapfile, "r"); - if (fp == NULL){ - sprintf(line,"Cannot open input map file %s", mapfile); - error->all(FLERR,line); - } - - if (fgets(line,MAXLINE,fp) == NULL) error->all(FLERR,"Error while reading header of mapping file!"); - nx = force->inumeric(FLERR, strtok(line, " \n\t\r\f")); - ny = force->inumeric(FLERR, strtok(NULL, " \n\t\r\f")); - nz = force->inumeric(FLERR, strtok(NULL, " \n\t\r\f")); - nucell = force->inumeric(FLERR, strtok(NULL, " \n\t\r\f")); - ntotal = nx*ny*nz; - if (ntotal*nucell != ngroup) error->all(FLERR,"FFT mesh and number of atoms in group mismatch!"); - - // second line of mapfile is comment - if (fgets(line,MAXLINE,fp) == NULL) error->all(FLERR,"Error while reading comment of mapping file!"); - - int ix, iy, iz, iu; - // the remaining lines carry the mapping info - for (int i = 0; i < ngroup; ++i){ - if (fgets(line,MAXLINE,fp) == NULL) {info = 1; break;} - ix = force->inumeric(FLERR, strtok(line, " \n\t\r\f")); - iy = force->inumeric(FLERR, strtok(NULL, " \n\t\r\f")); - iz = force->inumeric(FLERR, strtok(NULL, " \n\t\r\f")); - iu = force->inumeric(FLERR, strtok(NULL, " \n\t\r\f")); - itag = force->inumeric(FLERR, strtok(NULL, " \n\t\r\f")); - - // check if index is in correct range - if (ix < 0 || ix >= nx || iy < 0 || iy >= ny || iz < 0 || iz >= nz|| iu < 0 || iu >= nucell) {info = 2; break;} - // 1 <= itag <= natoms - if (itag < 1 || itag > static_cast(atom->natoms)) {info = 3; break;} - idx = ((ix*ny+iy)*nz+iz)*nucell + iu; - tag2surf[itag] = idx; - surf2tag[idx] = itag; - } - fclose(fp); - - if (tag2surf.size() != surf2tag.size() || tag2surf.size() != static_cast(ngroup) ) - error->all(FLERR,"The mapping is incomplete!"); - if (info) error->all(FLERR,"Error while reading mapping file!"); - - // check the correctness of mapping - int *mask = atom->mask; - int *tag = atom->tag; - int nlocal = atom->nlocal; - - for (int i = 0; i < nlocal; ++i) { - if (mask[i] & groupbit){ - itag = tag[i]; - idx = tag2surf[itag]; - if (itag != surf2tag[idx]) error->one(FLERR,"The mapping info read is incorrect!"); - } - } - -return; -} - -/* ---------------------------------------------------------------------- - * private method, to output the force constant matrix - * --------------------------------------------------------------------*/ -void FixPhonon::postprocess( ) -{ - if (neval < 1) return; - - ifreq = 0; - int idim, jdim, ndim; - double inv_neval = 1.0 /double(neval); - - // to get - for (idq = 0; idq < mynq; ++idq) - for (idim = 0; idim < fft_dim2; ++idim) Phi_q[idq][idim] = Rqsum[idq][idim]*inv_neval; - - // to get - for (idx = 0; idx < mynpt; ++idx) - for (idim = 0; idim < fft_dim; ++idim) Rnow[idx][idim] = Rsum[idx][idim] * inv_neval; - - // to get q - for (idim = 0; idim < fft_dim; ++idim){ - int m = 0; - for (idx = 0; idx < mynpt; ++idx){ - fft_data[m++] = Rnow[idx][idim]; - fft_data[m++] = 0.; - } - - fft->compute(fft_data,fft_data,-1); - - m = 0; - for (idq = 0; idq < mynq; ++idq){ - Rqnow[idq][idim] = std::complex(fft_data[m], fft_data[m+1]); - m += 2; - } - } - - // to get G(q) = - q.q - for (idq = 0; idq < mynq; ++idq){ - ndim = 0; - for (idim = 0; idim < fft_dim; ++idim) - for (jdim = 0; jdim < fft_dim; ++jdim) Phi_q[idq][ndim++] -= Rqnow[idq][idim]*conj(Rqnow[idq][jdim]); - } - - // to get Phi = KT.G^-1; normalization of FFTW data is done here - double boltz = force->boltz, kbtsqrt[sysdim], TempAve = 0.; - double TempFac = inv_neval*inv_nTemp; - double NormFac = TempFac*double(ntotal); - - for (idim = 0; idim < sysdim; ++idim){ - kbtsqrt[idim] = sqrt(TempSum[idim]*NormFac); - TempAve += TempSum[idim]*TempFac; - } - TempAve /= sysdim*boltz; - - for (idq = 0; idq < mynq; ++idq){ - GaussJordan(fft_dim, Phi_q[idq]); - ndim =0; - for (idim = 0; idim < fft_dim; ++idim) - for (jdim = 0; jdim < fft_dim; ++jdim) Phi_q[idq][ndim++] *= kbtsqrt[idim%sysdim]*kbtsqrt[jdim%sysdim]; - } - - // to collect all local Phi_q to root - displs[0]=0; - for (int i = 0; i < nprocs; ++i) recvcnts[i] = fft_cnts[i]*fft_dim*2; - for (int i = 1; i < nprocs; ++i) displs[i] = displs[i-1] + recvcnts[i-1]; - MPI_Gatherv(Phi_q[0],mynq*fft_dim2*2,MPI_DOUBLE,Phi_all[0],recvcnts,displs,MPI_DOUBLE,0,world); - - // to collect all basis info and averaged it on root - double basis_root[fft_dim]; - if (fft_dim > sysdim) MPI_Reduce(&basis[1][0], &basis_root[sysdim], fft_dim-sysdim, MPI_DOUBLE, MPI_SUM, 0, world); - - if (me == 0){ // output dynamic matrix by root - - // get basis info - for (idim = 0; idim < sysdim; ++idim) basis_root[idim] = 0.; - for (idim = sysdim; idim < fft_dim; ++idim) basis_root[idim] /= double(ntotal)*double(neval); - // get unit cell base vector info; might be incorrect if MD pbc and FixPhonon pbc mismatch. - double basevec[9]; - basevec[1] = basevec[2] = basevec[5] = 0.; - basevec[0] = hsum[0] * inv_neval / double(nx); - basevec[4] = hsum[1] * inv_neval / double(ny); - basevec[8] = hsum[2] * inv_neval / double(nz); - basevec[7] = hsum[3] * inv_neval / double(nz); - basevec[6] = hsum[4] * inv_neval / double(nz); - basevec[3] = hsum[5] * inv_neval / double(ny); - - // write binary file, in fact, it is the force constants matrix that is written - // Enforcement of ASR and the conversion of dynamical matrix is done in the postprocessing code - char fname[MAXLINE]; - sprintf(fname,"%s.bin." BIGINT_FORMAT,prefix,update->ntimestep); - FILE *fp_bin = fopen(fname,"wb"); - - fwrite(&sysdim, sizeof(int), 1, fp_bin); - fwrite(&nx, sizeof(int), 1, fp_bin); - fwrite(&ny, sizeof(int), 1, fp_bin); - fwrite(&nz, sizeof(int), 1, fp_bin); - fwrite(&nucell, sizeof(int), 1, fp_bin); - fwrite(&boltz, sizeof(double), 1, fp_bin); - - fwrite(Phi_all[0],sizeof(double),ntotal*fft_dim2*2,fp_bin); - - fwrite(&TempAve, sizeof(double),1, fp_bin); - fwrite(&basevec[0], sizeof(double),9, fp_bin); - fwrite(&basis_root[0],sizeof(double),fft_dim,fp_bin); - fwrite(basetype, sizeof(int), nucell, fp_bin); - fwrite(M_inv_sqrt, sizeof(double),nucell, fp_bin); - - fclose(fp_bin); - - // write log file, here however, it is the dynamical matrix that is written - for (int i = 0; i < 60; ++i) fprintf(flog,"#"); fprintf(flog,"\n"); - fprintf(flog, "# Current time step : " BIGINT_FORMAT "\n", update->ntimestep); - fprintf(flog, "# Total number of measurements : %d\n", neval); - fprintf(flog, "# Average temperature of the measurement : %lg\n", TempAve); - fprintf(flog, "# Boltzmann constant under current units : %lg\n", boltz); - fprintf(flog, "# basis vector A1 = [%lg %lg %lg]\n", basevec[0], basevec[1], basevec[2]); - fprintf(flog, "# basis vector A2 = [%lg %lg %lg]\n", basevec[3], basevec[4], basevec[5]); - fprintf(flog, "# basis vector A3 = [%lg %lg %lg]\n", basevec[6], basevec[7], basevec[8]); - for (int i = 0; i < 60; ++i) fprintf(flog,"#"); fprintf(flog,"\n"); - fprintf(flog, "# qx\t qy \t qz \t\t Phi(q)\n"); - - EnforceASR(); - - // to get D = 1/M x Phi - for (idq = 0; idq < ntotal; ++idq){ - ndim =0; - for (idim = 0; idim < fft_dim; ++idim) - for (jdim = 0; jdim < fft_dim; ++jdim) Phi_all[idq][ndim++] *= M_inv_sqrt[idim/sysdim]*M_inv_sqrt[jdim/sysdim]; - } - - idq =0; - for (int ix = 0; ix < nx; ++ix){ - double qx = double(ix)/double(nx); - for (int iy = 0; iy < ny; ++iy){ - double qy = double(iy)/double(ny); - for (int iz = 0; iz < nz; ++iz){ - double qz = double(iz)/double(nz); - fprintf(flog,"%lg %lg %lg", qx, qy, qz); - for (idim = 0; idim < fft_dim2; ++idim) fprintf(flog, " %lg %lg", real(Phi_all[idq][idim]), imag(Phi_all[idq][idim])); - fprintf(flog, "\n"); - ++idq; - } - } - } - fflush(flog); - } - -return; -} // end of postprocess - -/* ---------------------------------------------------------------------- - * private method, to get the inverse of a complex matrix by means of - * Gaussian-Jordan Elimination with full pivoting; square matrix required. - * - * Adapted from the Numerical Recipes in Fortran. - * --------------------------------------------------------------------*/ -void FixPhonon::GaussJordan(int n, std::complex *Mat) -{ - int i,icol,irow,j,k,l,ll,idr,idc; - int *indxc,*indxr,*ipiv; - double big, nmjk; - std::complex dum, pivinv; - - indxc = new int[n]; - indxr = new int[n]; - ipiv = new int[n]; - - for (i = 0; i < n; ++i) ipiv[i] = 0; - for (i = 0; i < n; ++i){ - big = 0.; - for (j = 0; j < n; ++j){ - if (ipiv[j] != 1){ - for (k = 0; k < n; ++k){ - if (ipiv[k] == 0){ - idr = j*n+k; - nmjk = norm(Mat[idr]); - if (nmjk >= big){ - big = nmjk; - irow = j; - icol = k; - } - } else if (ipiv[k] > 1) error->one(FLERR,"Singular matrix in complex GaussJordan!"); - } - } - } - ipiv[icol] += 1; - if (irow != icol){ - for (l = 0; l < n; ++l){ - idr = irow*n+l; - idc = icol*n+l; - dum = Mat[idr]; - Mat[idr] = Mat[idc]; - Mat[idc] = dum; - } - } - indxr[i] = irow; - indxc[i] = icol; - idr = icol*n+icol; - if (Mat[idr] == std::complex(0.,0.)) error->one(FLERR,"Singular matrix in complex GaussJordan!"); - - pivinv = 1./ Mat[idr]; - Mat[idr] = std::complex(1.,0.); - idr = icol*n; - for (l = 0; l < n; ++l) Mat[idr+l] *= pivinv; - for (ll = 0; ll < n; ++ll){ - if (ll != icol){ - idc = ll*n + icol; - dum = Mat[idc]; - Mat[idc] = 0.; - idc -= icol; - for (l = 0; l < n; ++l) Mat[idc+l] -= Mat[idr+l]*dum; - } - } - } - - for (l = n-1; l >= 0; --l){ - int rl = indxr[l]; - int cl = indxc[l]; - if (rl != cl){ - for (k = 0; k < n; ++k){ - idr = k*n + rl; - idc = k*n + cl; - dum = Mat[idr]; - Mat[idr] = Mat[idc]; - Mat[idc] = dum; - } - } - } - delete []indxr; - delete []indxc; - delete []ipiv; - -return; -} - -/* ---------------------------------------------------------------------- - * private method, to apply the acoustic sum rule on force constant matrix - * at gamma point. Should be executed on root only. - * --------------------------------------------------------------------*/ -void FixPhonon::EnforceASR() -{ - if (nasr < 1) return; - - for (int iit = 0; iit < nasr; ++iit){ - // simple ASR; the resultant matrix might not be symmetric - for (int a = 0; a < sysdim; ++a) - for (int b = 0; b < sysdim; ++b){ - for (int k = 0; k < nucell; ++k){ - double sum = 0.; - for (int kp = 0; kp < nucell; ++kp){ - int idx = (k*sysdim+a)*fft_dim + kp*sysdim + b; - sum += real(Phi_all[0][idx]); - } - sum /= double(nucell); - for (int kp = 0; kp < nucell; ++kp){ - int idx = (k*sysdim+a)*fft_dim + kp*sysdim + b; - real(Phi_all[0][idx]) -= sum; - } - } - } - - // symmetrize - for (int k = 0; k < nucell; ++k) - for (int kp = k; kp < nucell; ++kp){ - double csum = 0.; - for (int a = 0; a < sysdim; ++a) - for (int b = 0; b < sysdim; ++b){ - int idx = (k*sysdim+a)*fft_dim + kp*sysdim + b; - int jdx = (kp*sysdim+b)*fft_dim + k*sysdim + a; - csum = (real(Phi_all[0][idx])+real(Phi_all[0][jdx]))*0.5; - real(Phi_all[0][idx]) = real(Phi_all[0][jdx]) = csum; - } - } - } - - // symmetric ASR - for (int a = 0; a < sysdim; ++a) - for (int b = 0; b < sysdim; ++b){ - for (int k = 0; k < nucell; ++k){ - double sum = 0.; - for (int kp = 0; kp < nucell; ++kp){ - int idx = (k*sysdim+a)*fft_dim + kp*sysdim + b; - sum += real(Phi_all[0][idx]); - } - sum /= double(nucell-k); - for (int kp = k; kp < nucell; ++kp){ - int idx = (k*sysdim+a)*fft_dim + kp*sysdim + b; - int jdx = (kp*sysdim+b)*fft_dim + k*sysdim + a; - real(Phi_all[0][idx]) -= sum; - real(Phi_all[0][jdx]) = real(Phi_all[0][idx]); - } - } - } - -return; -} -/* --------------------------------------------------------------------*/ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: + Ling-Ti Kong + + Contact: + School of Materials Science and Engineering, + Shanghai Jiao Tong University, + 800 Dongchuan Road, Minhang, + Shanghai 200240, CHINA + + konglt@sjtu.edu.cn; konglt@gmail.com +------------------------------------------------------------------------- */ + +#include "string.h" +#include "fix_phonon.h" +#include "atom.h" +#include "compute.h" +#include "domain.h" +#include "fft3d_wrap.h" +#include "force.h" +#include "group.h" +#include "lattice.h" +#include "modify.h" +#include "update.h" +#include "citeme.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +#define INVOKED_SCALAR 1 +#define INVOKED_VECTOR 2 +#define MAXLINE 512 + +static const char cite_fix_phonon[] = + "fix phonon command:\n\n" + "@Article{Kong11,\n" + " author = {L. T. Kong},\n" + " title = {Phonon dispersion measured directly from molecular dynamics simulations},\n" + " journal = {Comp.~Phys.~Comm.},\n" + " year = 2011,\n" + " volume = 182,\n" + " pages = {2201--2207}\n" + "}\n\n"; + +/* ---------------------------------------------------------------------- */ + +FixPhonon::FixPhonon(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) +{ + if (lmp->citeme) lmp->citeme->add(cite_fix_phonon); + + MPI_Comm_rank(world,&me); + MPI_Comm_size(world,&nprocs); + + if (narg < 8) error->all(FLERR,"Illegal fix phonon command: number of arguments < 8"); + + nevery = force->inumeric(FLERR, arg[3]); // Calculate this fix every n steps! + if (nevery < 1) error->all(FLERR,"Illegal fix phonon command"); + + nfreq = force->inumeric(FLERR, arg[4]); // frequency to output result + if (nfreq < 1) error->all(FLERR,"Illegal fix phonon command"); + + waitsteps = ATOBIGINT(arg[5]); // Wait this many timesteps before actually measuring + if (waitsteps < 0) error->all(FLERR,"Illegal fix phonon command: waitsteps < 0 !"); + + int n = strlen(arg[6]) + 1; // map file + mapfile = new char[n]; + strcpy(mapfile, arg[6]); + + n = strlen(arg[7]) + 1; // prefix of output + prefix = new char[n]; + strcpy(prefix, arg[7]); + logfile = new char[n+4]; + sprintf(logfile,"%s.log",prefix); + + int sdim = sysdim = domain->dimension; + int iarg = 8; + nasr = 20; + + // other command line options + while (iarg < narg){ + if (strcmp(arg[iarg],"sysdim") == 0){ + if (++iarg >= narg) error->all(FLERR,"Illegal fix phonon command: incomplete command line options."); + sdim = force->inumeric(FLERR, arg[iarg]); + if (sdim < 1) error->all(FLERR,"Illegal fix phonon command: sysdim should not be less than 1."); + + } else if (strcmp(arg[iarg],"nasr") == 0){ + if (++iarg >= narg) error->all(FLERR,"Illegal fix phonon command: incomplete command line options."); + nasr = force->inumeric(FLERR, arg[iarg]); + + } else { + error->all(FLERR,"Illegal fix phonon command: unknown option read!"); + } + + ++iarg; + } + + // get the dimension of the simulation; 1D is possible by specifying the option of "sysdim 1" + if (sdim < sysdim) sysdim = sdim; + nasr = MAX(0, nasr); + + // get the total number of atoms in group + ngroup = static_cast(group->count(igroup)); + if (ngroup < 1) error->all(FLERR,"No atom found for fix phonon!"); + + // MPI gatherv related variables + recvcnts = new int[nprocs]; + displs = new int[nprocs]; + + // mapping index + tag2surf.clear(); // clear map info + surf2tag.clear(); + + // get the mapping between lattice indices and atom IDs + readmap(); delete []mapfile; + if (nucell == 1) nasr = MIN(1,nasr); + + // get the mass matrix for dynamic matrix + getmass(); + + // create FFT and allocate memory for FFT + // here the parallization is done on the x direction only + nxlo = 0; + int *nx_loc = new int [nprocs]; + for (int i = 0; i < nprocs; ++i){ + nx_loc[i] = nx/nprocs; + if (i < nx%nprocs) ++nx_loc[i]; + } + for (int i = 0; i < me; ++i) nxlo += nx_loc[i]; + nxhi = nxlo + nx_loc[me] - 1; + mynpt = nx_loc[me]*ny*nz; + mynq = mynpt; + + fft_dim = nucell*sysdim; + fft_dim2 = fft_dim*fft_dim; + fft_nsend = mynpt*fft_dim; + + fft_cnts = new int[nprocs]; + fft_disp = new int[nprocs]; + fft_disp[0] = 0; + for (int i = 0; i < nprocs; ++i) fft_cnts[i] = nx_loc[i]*ny*nz*fft_dim; + for (int i = 1; i < nprocs; ++i) fft_disp[i] = fft_disp[i-1]+fft_cnts[i-1]; + delete []nx_loc; + + fft = new FFT3d(lmp,world,nz,ny,nx,0,nz-1,0,ny-1,nxlo,nxhi,0,nz-1,0,ny-1,nxlo,nxhi,0,0,&mysize); + memory->create(fft_data, MAX(1,mynq)*2, "fix_phonon:fft_data"); + + // allocate variables; MAX(1,... is used because NULL buffer will result in error for MPI + memory->create(RIloc,ngroup,(sysdim+1),"fix_phonon:RIloc"); + memory->create(RIall,ngroup,(sysdim+1),"fix_phonon:RIall"); + memory->create(Rsort,ngroup, sysdim, "fix_phonon:Rsort"); + + memory->create(Rnow, MAX(1,mynpt),fft_dim,"fix_phonon:Rnow"); + memory->create(Rsum, MAX(1,mynpt),fft_dim,"fix_phonon:Rsum"); + + memory->create(basis,nucell, sysdim, "fix_phonon:basis"); + + // because of hermit, only nearly half of q points are stored + memory->create(Rqnow,MAX(1,mynq),fft_dim, "fix_phonon:Rqnow"); + memory->create(Rqsum,MAX(1,mynq),fft_dim2,"fix_phonon:Rqsum"); + memory->create(Phi_q,MAX(1,mynq),fft_dim2,"fix_phonon:Phi_q"); + + // variable to collect all local Phi to root + if (me == 0) memory->create(Phi_all,ntotal,fft_dim2,"fix_phonon:Phi_all"); + else memory->create(Phi_all,1,1,"fix_phonon:Phi_all"); + + // output some information on the system to log file + if (me == 0){ + flog = fopen(logfile, "w"); + if (flog == NULL) { + char str[MAXLINE]; + sprintf(str,"Can not open output file %s",logfile); + error->one(FLERR,str); + } + for (int i = 0; i < 60; ++i) fprintf(flog,"#"); fprintf(flog,"\n"); + fprintf(flog,"# group name of the atoms under study : %s\n", group->names[igroup]); + fprintf(flog,"# total number of atoms in the group : %d\n", ngroup); + fprintf(flog,"# dimension of the system : %d D\n", sysdim); + fprintf(flog,"# number of atoms per unit cell : %d\n", nucell); + fprintf(flog,"# dimension of the FFT mesh : %d x %d x %d\n", nx, ny, nz); + fprintf(flog,"# number of wait steps before measurement : " BIGINT_FORMAT "\n", waitsteps); + fprintf(flog,"# frequency of the measurement : %d\n", nevery); + fprintf(flog,"# output result after this many measurement: %d\n", nfreq); + fprintf(flog,"# number of processors used by this run : %d\n", nprocs); + for (int i = 0; i < 60; ++i) fprintf(flog,"#"); fprintf(flog,"\n"); + fprintf(flog,"# mapping information between lattice index and atom id\n"); + fprintf(flog,"# nx ny nz nucell\n"); + fprintf(flog,"%d %d %d %d\n", nx, ny, nz, nucell); + fprintf(flog,"# l1 l2 l3 k atom_id\n"); + int ix, iy, iz, iu; + for (idx = 0; idx < ngroup; ++idx){ + itag = surf2tag[idx]; + iu = idx%nucell; + iz = (idx/nucell)%nz; + iy = (idx/(nucell*nz))%ny; + ix = (idx/(nucell*nz*ny))%nx; + fprintf(flog,"%d %d %d %d %d\n", ix, iy, iz, iu, itag); + } + for (int i = 0; i < 60; ++i) fprintf(flog,"#"); fprintf(flog,"\n"); + fflush(flog); + } + surf2tag.clear(); + + // default temperature is from thermo + TempSum = new double[sysdim]; + id_temp = new char[12]; + strcpy(id_temp,"thermo_temp"); + int icompute = modify->find_compute(id_temp); + temperature = modify->compute[icompute]; + inv_nTemp = 1.0/group->count(temperature->igroup); + +return; +} // end of constructor + +/* ---------------------------------------------------------------------- */ + +void FixPhonon::post_run() +{ + // compute and output final results + if (ifreq > 0 && ifreq != nfreq) postprocess(); + if (me == 0) fclose(flog); +} + +/* ---------------------------------------------------------------------- */ + +FixPhonon::~FixPhonon() +{ + // delete locally stored array + memory->destroy(RIloc); + memory->destroy(RIall); + memory->destroy(Rsort); + memory->destroy(Rnow); + memory->destroy(Rsum); + + memory->destroy(basis); + + memory->destroy(Rqnow); + memory->destroy(Rqsum); + memory->destroy(Phi_q); + memory->destroy(Phi_all); + + delete []recvcnts; + delete []displs; + delete []prefix; + delete []logfile; + delete []fft_cnts; + delete []fft_disp; + delete []id_temp; + delete []TempSum; + delete []M_inv_sqrt; + delete []basetype; + + // destroy FFT + delete fft; + memory->sfree(fft_data); + + // clear map info + tag2surf.clear(); + surf2tag.clear(); + +} + +/* ---------------------------------------------------------------------- */ + +int FixPhonon::setmask() +{ + int mask = 0; + mask |= END_OF_STEP; + +return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixPhonon::init() +{ + // warn if more than one fix-phonon + int count = 0; + for (int i = 0; i < modify->nfix; ++i) if (strcmp(modify->fix[i]->style,"gfc") == 0) ++count; + if (count > 1 && me == 0) error->warning(FLERR,"More than one fix phonon defined"); // just warn, but allowed. + +return; +} + +/* ---------------------------------------------------------------------- */ + +void FixPhonon::setup(int flag) +{ + // initialize accumulating variables + for (int i = 0; i < sysdim; ++i) TempSum[i] = 0.; + + for (int i = 0; i < mynpt; ++i) + for (int j = 0; j < fft_dim; ++j) Rsum[i][j] = 0.; + + for (int i =0; i < mynq; ++i) + for (int j =0; j < fft_dim2; ++j) Rqsum[i][j] = std::complex (0.,0.); + + for (int i = 0; i < 6; ++i) hsum[i] = 0.; + + for (int i = 0; i < nucell; ++i) + for (int j = 0; j < sysdim; ++j) basis[i][j] = 0.; + + prev_nstep = update->ntimestep; + neval = ifreq = 0; + +return; +} + +/* ---------------------------------------------------------------------- */ + +void FixPhonon::end_of_step() +{ + if ( (update->ntimestep-prev_nstep) <= waitsteps) return; + + double **x = atom->x; + int *mask = atom->mask; + int *tag = atom->tag; + int *image = atom->image; + int nlocal = atom->nlocal; + + double xprd = domain->xprd; + double yprd = domain->yprd; + double zprd = domain->zprd; + double *h = domain->h; + double xbox, ybox, zbox; + + int i,idim,jdim,ndim; + double xcur[3]; + + // to get the current temperature + if (!(temperature->invoked_flag & INVOKED_VECTOR)) temperature->compute_vector(); + for (idim = 0; idim < sysdim; ++idim) TempSum[idim] += temperature->vector[idim]; + + // evaluate R(r) on local proc + nfind = 0; + for (i = 0; i < nlocal; ++i){ + if (mask[i] & groupbit){ + itag = tag[i]; + idx = tag2surf[itag]; + + domain->unmap(x[i], image[i], xcur); + + for (idim = 0; idim < sysdim; ++idim) RIloc[nfind][idim] = xcur[idim]; + RIloc[nfind++][sysdim] = double(idx); + } + } + + // gather R(r) on local proc, then sort and redistribute to all procs for FFT + nfind *= (sysdim+1); + displs[0] = 0; + for (i = 0; i < nprocs; ++i) recvcnts[i] = 0; + MPI_Gather(&nfind,1,MPI_INT,recvcnts,1,MPI_INT,0,world); + for (i = 1; i < nprocs; ++i) displs[i] = displs[i-1] + recvcnts[i-1]; + + MPI_Gatherv(RIloc[0],nfind,MPI_DOUBLE,RIall[0],recvcnts,displs,MPI_DOUBLE,0,world); + if (me == 0){ + for (i = 0; i < ngroup; ++i){ + idx = static_cast(RIall[i][sysdim]); + for (idim = 0; idim < sysdim; ++idim) Rsort[idx][idim] = RIall[i][idim]; + } + } + MPI_Scatterv(Rsort[0],fft_cnts,fft_disp, MPI_DOUBLE, Rnow[0], fft_nsend, MPI_DOUBLE,0,world); + + // get Rsum + for (idx = 0; idx < mynpt; ++idx) + for (idim = 0; idim < fft_dim; ++idim) Rsum[idx][idim] += Rnow[idx][idim]; + + // FFT R(r) to get R(q) + for (idim = 0; idim < fft_dim; ++idim){ + int m=0; + for (idx = 0; idx < mynpt; ++idx){ + fft_data[m++] = Rnow[idx][idim]; + fft_data[m++] = 0.; + } + + fft->compute(fft_data, fft_data, -1); + + m = 0; + for (idq = 0; idq < mynq; ++idq){ + Rqnow[idq][idim] = std::complex(fft_data[m], fft_data[m+1]); + m += 2; + } + } + + // to get sum(R(q).R(q)*) + for (idq = 0; idq < mynq; ++idq){ + ndim = 0; + for (idim = 0; idim < fft_dim; ++idim) + for (jdim = 0; jdim < fft_dim; ++jdim) Rqsum[idq][ndim++] += Rqnow[idq][idim]*conj(Rqnow[idq][jdim]); + } + + // get basis info + if (fft_dim > sysdim){ + double dist2orig[3]; + for (idx = 0; idx < mynpt; ++idx){ + ndim = sysdim; + for (i = 1; i < nucell; ++i){ + for (idim = 0; idim < sysdim; ++idim) dist2orig[idim] = Rnow[idx][ndim++] - Rnow[idx][idim]; + domain->minimum_image(dist2orig); + for (idim = 0; idim < sysdim; ++idim) basis[i][idim] += dist2orig[idim]; + } + } + } + // get lattice vector info + for (int i = 0; i < 6; ++i) hsum[i] += h[i]; + + // increment counter + ++neval; + + // compute and output Phi_q after every nfreq evaluations + if (++ifreq == nfreq) postprocess(); + +return; +} // end of end_of_step() + +/* ---------------------------------------------------------------------- */ + +double FixPhonon::memory_usage() +{ + double bytes = sizeof(double)*2*mynq + + sizeof(std::map)*2*ngroup + + sizeof(double)*(ngroup*(3*sysdim+2)+mynpt*fft_dim*2) + + sizeof(std::complex)*MAX(1,mynq)*fft_dim *(1+2*fft_dim) + + sizeof(std::complex)*ntotal*fft_dim2 + + sizeof(int) * nprocs * 4; + return bytes; +} + +/* ---------------------------------------------------------------------- */ + +int FixPhonon::modify_param(int narg, char **arg) +{ + if (strcmp(arg[0],"temp") == 0) { + if (narg < 2) error->all(FLERR,"Illegal fix_modify command"); + delete [] id_temp; + int n = strlen(arg[1]) + 1; + id_temp = new char[n]; + strcpy(id_temp,arg[1]); + + int icompute = modify->find_compute(id_temp); + if (icompute < 0) error->all(FLERR,"Could not find fix_modify temp ID"); + temperature = modify->compute[icompute]; + + if (temperature->tempflag == 0) + error->all(FLERR,"Fix_modify temp ID does not compute temperature"); + inv_nTemp = 1.0/group->count(temperature->igroup); + + return 2; + } + +return 0; +} + +/* ---------------------------------------------------------------------- + * private method, to get the mass matrix for dynamic matrix + * --------------------------------------------------------------------*/ +void FixPhonon::getmass() +{ + int nlocal = atom->nlocal; + int *mask = atom->mask; + int *tag = atom->tag; + int *type = atom->type; + double *rmass = atom->rmass; + double *mass = atom->mass; + double *mass_one, *mass_all; + double *type_one, *type_all; + + mass_one = new double[nucell]; + mass_all = new double[nucell]; + type_one = new double[nucell]; + type_all = new double[nucell]; + for (int i = 0; i < nucell; ++i) mass_one[i] = type_one[i] = 0.; + + if (rmass){ + for (int i = 0; i < nlocal; ++i){ + if (mask[i] & groupbit){ + itag = tag[i]; + idx = tag2surf[itag]; + int iu = idx%nucell; + mass_one[iu] += rmass[i]; + type_one[iu] += double(type[i]); + } + } + } else { + for (int i = 0; i < nlocal; ++i){ + if (mask[i] & groupbit){ + itag = tag[i]; + idx = tag2surf[itag]; + int iu = idx%nucell; + mass_one[iu] += mass[type[i]]; + type_one[iu] += double(type[i]); + } + } + } + + MPI_Allreduce(mass_one,mass_all,nucell,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(type_one,type_all,nucell,MPI_DOUBLE,MPI_SUM,world); + + M_inv_sqrt = new double[nucell]; + basetype = new int[nucell]; + + double inv_total = 1./double(ntotal); + for (int i = 0; i < nucell; ++i){ + mass_all[i] *= inv_total; + M_inv_sqrt[i] = sqrt(1./mass_all[i]); + + basetype[i] = int(type_all[i]*inv_total); + } + delete []mass_one; + delete []mass_all; + delete []type_one; + delete []type_all; + +return; +} + + +/* ---------------------------------------------------------------------- + * private method, to read the mapping info from file + * --------------------------------------------------------------------*/ + +void FixPhonon::readmap() +{ + int info = 0; + + // auto-generate mapfile for "cluster" (gamma only system) + if (strcmp(mapfile, "GAMMA") == 0){ + nx = ny = nz = ntotal = 1; + nucell = ngroup; + + int *tag_loc, *tag_all; + memory->create(tag_loc,ngroup,"fix_phonon:tag_loc"); + memory->create(tag_all,ngroup,"fix_phonon:tag_all"); + + // get atom IDs on local proc + int nfind = 0; + for (int i = 0; i < atom->nlocal; ++i){ + if (atom->mask[i] & groupbit) tag_loc[nfind++] = atom->tag[i]; + } + + // gather IDs on local proc + displs[0] = 0; + for (int i = 0; i < nprocs; ++i) recvcnts[i] = 0; + MPI_Allgather(&nfind,1,MPI_INT,recvcnts,1,MPI_INT,world); + for (int i = 1; i < nprocs; ++i) displs[i] = displs[i-1] + recvcnts[i-1]; + + MPI_Allgatherv(tag_loc,nfind,MPI_INT,tag_all,recvcnts,displs,MPI_INT,world); + for (int i = 0; i < ngroup; ++i){ + itag = tag_all[i]; + tag2surf[itag] = i; + surf2tag[i] = itag; + } + + memory->destroy(tag_loc); + memory->destroy(tag_all); + return; + } + + // read from map file for others + char line[MAXLINE]; + FILE *fp = fopen(mapfile, "r"); + if (fp == NULL){ + sprintf(line,"Cannot open input map file %s", mapfile); + error->all(FLERR,line); + } + + if (fgets(line,MAXLINE,fp) == NULL) error->all(FLERR,"Error while reading header of mapping file!"); + nx = force->inumeric(FLERR, strtok(line, " \n\t\r\f")); + ny = force->inumeric(FLERR, strtok(NULL, " \n\t\r\f")); + nz = force->inumeric(FLERR, strtok(NULL, " \n\t\r\f")); + nucell = force->inumeric(FLERR, strtok(NULL, " \n\t\r\f")); + ntotal = nx*ny*nz; + if (ntotal*nucell != ngroup) error->all(FLERR,"FFT mesh and number of atoms in group mismatch!"); + + // second line of mapfile is comment + if (fgets(line,MAXLINE,fp) == NULL) error->all(FLERR,"Error while reading comment of mapping file!"); + + int ix, iy, iz, iu; + // the remaining lines carry the mapping info + for (int i = 0; i < ngroup; ++i){ + if (fgets(line,MAXLINE,fp) == NULL) {info = 1; break;} + ix = force->inumeric(FLERR, strtok(line, " \n\t\r\f")); + iy = force->inumeric(FLERR, strtok(NULL, " \n\t\r\f")); + iz = force->inumeric(FLERR, strtok(NULL, " \n\t\r\f")); + iu = force->inumeric(FLERR, strtok(NULL, " \n\t\r\f")); + itag = force->inumeric(FLERR, strtok(NULL, " \n\t\r\f")); + + // check if index is in correct range + if (ix < 0 || ix >= nx || iy < 0 || iy >= ny || iz < 0 || iz >= nz|| iu < 0 || iu >= nucell) {info = 2; break;} + // 1 <= itag <= natoms + if (itag < 1 || itag > static_cast(atom->natoms)) {info = 3; break;} + idx = ((ix*ny+iy)*nz+iz)*nucell + iu; + tag2surf[itag] = idx; + surf2tag[idx] = itag; + } + fclose(fp); + + if (tag2surf.size() != surf2tag.size() || tag2surf.size() != static_cast(ngroup) ) + error->all(FLERR,"The mapping is incomplete!"); + if (info) error->all(FLERR,"Error while reading mapping file!"); + + // check the correctness of mapping + int *mask = atom->mask; + int *tag = atom->tag; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; ++i) { + if (mask[i] & groupbit){ + itag = tag[i]; + idx = tag2surf[itag]; + if (itag != surf2tag[idx]) error->one(FLERR,"The mapping info read is incorrect!"); + } + } + +return; +} + +/* ---------------------------------------------------------------------- + * private method, to output the force constant matrix + * --------------------------------------------------------------------*/ +void FixPhonon::postprocess( ) +{ + if (neval < 1) return; + + ifreq = 0; + int idim, jdim, ndim; + double inv_neval = 1.0 /double(neval); + + // to get + for (idq = 0; idq < mynq; ++idq) + for (idim = 0; idim < fft_dim2; ++idim) Phi_q[idq][idim] = Rqsum[idq][idim]*inv_neval; + + // to get + for (idx = 0; idx < mynpt; ++idx) + for (idim = 0; idim < fft_dim; ++idim) Rnow[idx][idim] = Rsum[idx][idim] * inv_neval; + + // to get q + for (idim = 0; idim < fft_dim; ++idim){ + int m = 0; + for (idx = 0; idx < mynpt; ++idx){ + fft_data[m++] = Rnow[idx][idim]; + fft_data[m++] = 0.; + } + + fft->compute(fft_data,fft_data,-1); + + m = 0; + for (idq = 0; idq < mynq; ++idq){ + Rqnow[idq][idim] = std::complex(fft_data[m], fft_data[m+1]); + m += 2; + } + } + + // to get G(q) = - q.q + for (idq = 0; idq < mynq; ++idq){ + ndim = 0; + for (idim = 0; idim < fft_dim; ++idim) + for (jdim = 0; jdim < fft_dim; ++jdim) Phi_q[idq][ndim++] -= Rqnow[idq][idim]*conj(Rqnow[idq][jdim]); + } + + // to get Phi = KT.G^-1; normalization of FFTW data is done here + double boltz = force->boltz, kbtsqrt[sysdim], TempAve = 0.; + double TempFac = inv_neval*inv_nTemp; + double NormFac = TempFac*double(ntotal); + + for (idim = 0; idim < sysdim; ++idim){ + kbtsqrt[idim] = sqrt(TempSum[idim]*NormFac); + TempAve += TempSum[idim]*TempFac; + } + TempAve /= sysdim*boltz; + + for (idq = 0; idq < mynq; ++idq){ + GaussJordan(fft_dim, Phi_q[idq]); + ndim =0; + for (idim = 0; idim < fft_dim; ++idim) + for (jdim = 0; jdim < fft_dim; ++jdim) Phi_q[idq][ndim++] *= kbtsqrt[idim%sysdim]*kbtsqrt[jdim%sysdim]; + } + + // to collect all local Phi_q to root + displs[0]=0; + for (int i = 0; i < nprocs; ++i) recvcnts[i] = fft_cnts[i]*fft_dim*2; + for (int i = 1; i < nprocs; ++i) displs[i] = displs[i-1] + recvcnts[i-1]; + MPI_Gatherv(Phi_q[0],mynq*fft_dim2*2,MPI_DOUBLE,Phi_all[0],recvcnts,displs,MPI_DOUBLE,0,world); + + // to collect all basis info and averaged it on root + double basis_root[fft_dim]; + if (fft_dim > sysdim) MPI_Reduce(&basis[1][0], &basis_root[sysdim], fft_dim-sysdim, MPI_DOUBLE, MPI_SUM, 0, world); + + if (me == 0){ // output dynamic matrix by root + + // get basis info + for (idim = 0; idim < sysdim; ++idim) basis_root[idim] = 0.; + for (idim = sysdim; idim < fft_dim; ++idim) basis_root[idim] /= double(ntotal)*double(neval); + // get unit cell base vector info; might be incorrect if MD pbc and FixPhonon pbc mismatch. + double basevec[9]; + basevec[1] = basevec[2] = basevec[5] = 0.; + basevec[0] = hsum[0] * inv_neval / double(nx); + basevec[4] = hsum[1] * inv_neval / double(ny); + basevec[8] = hsum[2] * inv_neval / double(nz); + basevec[7] = hsum[3] * inv_neval / double(nz); + basevec[6] = hsum[4] * inv_neval / double(nz); + basevec[3] = hsum[5] * inv_neval / double(ny); + + // write binary file, in fact, it is the force constants matrix that is written + // Enforcement of ASR and the conversion of dynamical matrix is done in the postprocessing code + char fname[MAXLINE]; + sprintf(fname,"%s.bin." BIGINT_FORMAT,prefix,update->ntimestep); + FILE *fp_bin = fopen(fname,"wb"); + + fwrite(&sysdim, sizeof(int), 1, fp_bin); + fwrite(&nx, sizeof(int), 1, fp_bin); + fwrite(&ny, sizeof(int), 1, fp_bin); + fwrite(&nz, sizeof(int), 1, fp_bin); + fwrite(&nucell, sizeof(int), 1, fp_bin); + fwrite(&boltz, sizeof(double), 1, fp_bin); + + fwrite(Phi_all[0],sizeof(double),ntotal*fft_dim2*2,fp_bin); + + fwrite(&TempAve, sizeof(double),1, fp_bin); + fwrite(&basevec[0], sizeof(double),9, fp_bin); + fwrite(&basis_root[0],sizeof(double),fft_dim,fp_bin); + fwrite(basetype, sizeof(int), nucell, fp_bin); + fwrite(M_inv_sqrt, sizeof(double),nucell, fp_bin); + + fclose(fp_bin); + + // write log file, here however, it is the dynamical matrix that is written + for (int i = 0; i < 60; ++i) fprintf(flog,"#"); fprintf(flog,"\n"); + fprintf(flog, "# Current time step : " BIGINT_FORMAT "\n", update->ntimestep); + fprintf(flog, "# Total number of measurements : %d\n", neval); + fprintf(flog, "# Average temperature of the measurement : %lg\n", TempAve); + fprintf(flog, "# Boltzmann constant under current units : %lg\n", boltz); + fprintf(flog, "# basis vector A1 = [%lg %lg %lg]\n", basevec[0], basevec[1], basevec[2]); + fprintf(flog, "# basis vector A2 = [%lg %lg %lg]\n", basevec[3], basevec[4], basevec[5]); + fprintf(flog, "# basis vector A3 = [%lg %lg %lg]\n", basevec[6], basevec[7], basevec[8]); + for (int i = 0; i < 60; ++i) fprintf(flog,"#"); fprintf(flog,"\n"); + fprintf(flog, "# qx\t qy \t qz \t\t Phi(q)\n"); + + EnforceASR(); + + // to get D = 1/M x Phi + for (idq = 0; idq < ntotal; ++idq){ + ndim =0; + for (idim = 0; idim < fft_dim; ++idim) + for (jdim = 0; jdim < fft_dim; ++jdim) Phi_all[idq][ndim++] *= M_inv_sqrt[idim/sysdim]*M_inv_sqrt[jdim/sysdim]; + } + + idq =0; + for (int ix = 0; ix < nx; ++ix){ + double qx = double(ix)/double(nx); + for (int iy = 0; iy < ny; ++iy){ + double qy = double(iy)/double(ny); + for (int iz = 0; iz < nz; ++iz){ + double qz = double(iz)/double(nz); + fprintf(flog,"%lg %lg %lg", qx, qy, qz); + for (idim = 0; idim < fft_dim2; ++idim) fprintf(flog, " %lg %lg", real(Phi_all[idq][idim]), imag(Phi_all[idq][idim])); + fprintf(flog, "\n"); + ++idq; + } + } + } + fflush(flog); + } + +return; +} // end of postprocess + +/* ---------------------------------------------------------------------- + * private method, to get the inverse of a complex matrix by means of + * Gaussian-Jordan Elimination with full pivoting; square matrix required. + * + * Adapted from the Numerical Recipes in Fortran. + * --------------------------------------------------------------------*/ +void FixPhonon::GaussJordan(int n, std::complex *Mat) +{ + int i,icol,irow,j,k,l,ll,idr,idc; + int *indxc,*indxr,*ipiv; + double big, nmjk; + std::complex dum, pivinv; + + indxc = new int[n]; + indxr = new int[n]; + ipiv = new int[n]; + + for (i = 0; i < n; ++i) ipiv[i] = 0; + for (i = 0; i < n; ++i){ + big = 0.; + for (j = 0; j < n; ++j){ + if (ipiv[j] != 1){ + for (k = 0; k < n; ++k){ + if (ipiv[k] == 0){ + idr = j*n+k; + nmjk = norm(Mat[idr]); + if (nmjk >= big){ + big = nmjk; + irow = j; + icol = k; + } + } else if (ipiv[k] > 1) error->one(FLERR,"Singular matrix in complex GaussJordan!"); + } + } + } + ipiv[icol] += 1; + if (irow != icol){ + for (l = 0; l < n; ++l){ + idr = irow*n+l; + idc = icol*n+l; + dum = Mat[idr]; + Mat[idr] = Mat[idc]; + Mat[idc] = dum; + } + } + indxr[i] = irow; + indxc[i] = icol; + idr = icol*n+icol; + if (Mat[idr] == std::complex(0.,0.)) error->one(FLERR,"Singular matrix in complex GaussJordan!"); + + pivinv = 1./ Mat[idr]; + Mat[idr] = std::complex(1.,0.); + idr = icol*n; + for (l = 0; l < n; ++l) Mat[idr+l] *= pivinv; + for (ll = 0; ll < n; ++ll){ + if (ll != icol){ + idc = ll*n + icol; + dum = Mat[idc]; + Mat[idc] = 0.; + idc -= icol; + for (l = 0; l < n; ++l) Mat[idc+l] -= Mat[idr+l]*dum; + } + } + } + + for (l = n-1; l >= 0; --l){ + int rl = indxr[l]; + int cl = indxc[l]; + if (rl != cl){ + for (k = 0; k < n; ++k){ + idr = k*n + rl; + idc = k*n + cl; + dum = Mat[idr]; + Mat[idr] = Mat[idc]; + Mat[idc] = dum; + } + } + } + delete []indxr; + delete []indxc; + delete []ipiv; + +return; +} + +/* ---------------------------------------------------------------------- + * private method, to apply the acoustic sum rule on force constant matrix + * at gamma point. Should be executed on root only. + * --------------------------------------------------------------------*/ +void FixPhonon::EnforceASR() +{ + if (nasr < 1) return; + + for (int iit = 0; iit < nasr; ++iit){ + // simple ASR; the resultant matrix might not be symmetric + for (int a = 0; a < sysdim; ++a) + for (int b = 0; b < sysdim; ++b){ + for (int k = 0; k < nucell; ++k){ + double sum = 0.; + for (int kp = 0; kp < nucell; ++kp){ + int idx = (k*sysdim+a)*fft_dim + kp*sysdim + b; + sum += real(Phi_all[0][idx]); + } + sum /= double(nucell); + for (int kp = 0; kp < nucell; ++kp){ + int idx = (k*sysdim+a)*fft_dim + kp*sysdim + b; + real(Phi_all[0][idx]) -= sum; + } + } + } + + // symmetrize + for (int k = 0; k < nucell; ++k) + for (int kp = k; kp < nucell; ++kp){ + double csum = 0.; + for (int a = 0; a < sysdim; ++a) + for (int b = 0; b < sysdim; ++b){ + int idx = (k*sysdim+a)*fft_dim + kp*sysdim + b; + int jdx = (kp*sysdim+b)*fft_dim + k*sysdim + a; + csum = (real(Phi_all[0][idx])+real(Phi_all[0][jdx]))*0.5; + real(Phi_all[0][idx]) = real(Phi_all[0][jdx]) = csum; + } + } + } + + // symmetric ASR + for (int a = 0; a < sysdim; ++a) + for (int b = 0; b < sysdim; ++b){ + for (int k = 0; k < nucell; ++k){ + double sum = 0.; + for (int kp = 0; kp < nucell; ++kp){ + int idx = (k*sysdim+a)*fft_dim + kp*sysdim + b; + sum += real(Phi_all[0][idx]); + } + sum /= double(nucell-k); + for (int kp = k; kp < nucell; ++kp){ + int idx = (k*sysdim+a)*fft_dim + kp*sysdim + b; + int jdx = (kp*sysdim+b)*fft_dim + k*sysdim + a; + real(Phi_all[0][idx]) -= sum; + real(Phi_all[0][jdx]) = real(Phi_all[0][idx]); + } + } + } + +return; +} +/* --------------------------------------------------------------------*/ diff --git a/src/USER-PHONON/fix_phonon.h b/src/USER-PHONON/fix_phonon.h index 96bc04b91f..e9e0aa79ad 100644 --- a/src/USER-PHONON/fix_phonon.h +++ b/src/USER-PHONON/fix_phonon.h @@ -1,185 +1,185 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - www.cs.sandia.gov/~sjplimp/lammps.html - Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ -/* ---------------------------------------------------------------------- - Contributing authors: - Ling-Ti Kong - - Contact: - School of Materials Science and Engineering, - Shanghai Jiao Tong University, - 800 Dongchuan Road, Minhang, - Shanghai 200240, CHINA - - konglt@sjtu.edu.cn; konglt@gmail.com -------------------------------------------------------------------------- */ -#ifdef FIX_CLASS - -FixStyle(phonon,FixPhonon) - -#else - -#ifndef FIX_PHONON_H -#define FIX_PHONON_H - -#include -#include "fix.h" -#include -#include "stdio.h" -#include "stdlib.h" -#include "string.h" - -namespace LAMMPS_NS { - -class FixPhonon : public Fix { - public: - FixPhonon(class LAMMPS *, int, char **); - ~FixPhonon(); - - int setmask(); - void init(); - void setup(int); - void end_of_step(); - void post_run(); - double memory_usage(); - int modify_param(int, char **); - - private: - int me,nprocs; - bigint waitsteps; // wait these number of timesteps before recording atom positions - bigint prev_nstep; // number of steps from previous run(s); to judge if waitsteps is reached. - int nfreq, ifreq; // after this number of measurement (nfreq), the result will be output once - int nx,ny,nz,nucell,ntotal; // surface dimensions in x- and y-direction, number of atom per unit surface cell - int neval; // # of evaluations - int sysdim; // system dimension - int ngroup, nfind; // total number of atoms in group; total number of atoms on this proc - char *prefix, *logfile; // prefix of output file names - FILE *flog; - - double *M_inv_sqrt; - - class FFT3d *fft; // to do fft via the fft3d wraper - int nxlo,nxhi,mysize; // size info for local MPI_FFTW - int mynpt,mynq,fft_nsend; - int *fft_cnts, *fft_disp; - int fft_dim, fft_dim2; - double *fft_data; - - int itag, idx, idq; // index variables - std::map tag2surf, surf2tag; // Mapping info - - double **RIloc; // R(r) and index on local proc - double **RIall; // gathered R(r) and index - double **Rsort; // sorted R(r) - double **Rnow; // Current R(r) on local proc - double **Rsum; // Accumulated R(r) on local proc - - int *recvcnts, *displs; // MPI related variables - - std::complex **Rqnow; // Current R(q) on local proc - std::complex **Rqsum; // Accumulator for conj(R(q)_alpha)*R(q)_beta - std::complex **Phi_q; // Phi's on local proc - std::complex **Phi_all; // Phi for all - - void readmap(); // to read the mapping of gf atoms - char *mapfile; // file name of the map file - - void getmass(); // to get the mass of each atom in a unit cell - - int nasr; - void postprocess(); // to post process the data - void EnforceASR(); // to apply acoustic sum rule to gamma point force constant matrix - - char *id_temp; // compute id for temperature - double *TempSum; // to get the average temperature vector - double inv_nTemp; // inverse of number of atoms in temperature group - class Compute *temperature; // compute that computes the temperature - - double hsum[6], **basis; - int *basetype; - - // private methods to do matrix inversion - void GaussJordan(int, std::complex*); - -}; -} -#endif -#endif - -/* ERROR/WARNING messages: - -E: Illegal fix phonon command... - -Self-explanatory. Check the input script syntax and compare to the -documentation for the command. You can use -echo screen as a -command-line option when running LAMMPS to see the offending line. - -E: No atom found for fix phonon! - -Self-explanatory. Number of atoms in the group that was passed to -fix-phonon is less than 1. - -E: Can not open output file %s" - -Self-explanatory. - -E: Illegal fix_modify command - -Self-explanatory. - -E: Could not find fix_modify temp ID - -Self-explanatory. - -E: Fix_modify temp ID does not compute temperature - -Self-explanatory. - -E: Cannot open input map file %s - -Self-explanatory. - -E: Error while reading header of mapping file! - -Self-explanatory. The first line of the map file is expected to -contain 4 positive integer numbers. - -E: FFT mesh and number of atoms in group mismatch! - -Self-explanatory. The product of the 4 numbers should be exactly the -total number of atoms in the group that was passed to fix-phonon. - -E: Error while reading comment of mapping file! - -Self-explanatory. The second line of the map file should be a comment line. - -E: The mapping is incomplete! - -Self-explanatory. - -E: Error while reading mapping file! - -Self-explanatory. - -E: The mapping info read is incorrect! - -Self-explanatory. - -E: Singular matrix in complex GaussJordan! - -Self-explanatory. - -W: More than one fix phonon defined - -Self-explanatory. Just to warn that more than one fix-phonon is defined, but allowed. - -*/ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + www.cs.sandia.gov/~sjplimp/lammps.html + Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: + Ling-Ti Kong + + Contact: + School of Materials Science and Engineering, + Shanghai Jiao Tong University, + 800 Dongchuan Road, Minhang, + Shanghai 200240, CHINA + + konglt@sjtu.edu.cn; konglt@gmail.com +------------------------------------------------------------------------- */ +#ifdef FIX_CLASS + +FixStyle(phonon,FixPhonon) + +#else + +#ifndef FIX_PHONON_H +#define FIX_PHONON_H + +#include +#include "fix.h" +#include +#include "stdio.h" +#include "stdlib.h" +#include "string.h" + +namespace LAMMPS_NS { + +class FixPhonon : public Fix { + public: + FixPhonon(class LAMMPS *, int, char **); + ~FixPhonon(); + + int setmask(); + void init(); + void setup(int); + void end_of_step(); + void post_run(); + double memory_usage(); + int modify_param(int, char **); + + private: + int me,nprocs; + bigint waitsteps; // wait these number of timesteps before recording atom positions + bigint prev_nstep; // number of steps from previous run(s); to judge if waitsteps is reached. + int nfreq, ifreq; // after this number of measurement (nfreq), the result will be output once + int nx,ny,nz,nucell,ntotal; // surface dimensions in x- and y-direction, number of atom per unit surface cell + int neval; // # of evaluations + int sysdim; // system dimension + int ngroup, nfind; // total number of atoms in group; total number of atoms on this proc + char *prefix, *logfile; // prefix of output file names + FILE *flog; + + double *M_inv_sqrt; + + class FFT3d *fft; // to do fft via the fft3d wraper + int nxlo,nxhi,mysize; // size info for local MPI_FFTW + int mynpt,mynq,fft_nsend; + int *fft_cnts, *fft_disp; + int fft_dim, fft_dim2; + double *fft_data; + + int itag, idx, idq; // index variables + std::map tag2surf, surf2tag; // Mapping info + + double **RIloc; // R(r) and index on local proc + double **RIall; // gathered R(r) and index + double **Rsort; // sorted R(r) + double **Rnow; // Current R(r) on local proc + double **Rsum; // Accumulated R(r) on local proc + + int *recvcnts, *displs; // MPI related variables + + std::complex **Rqnow; // Current R(q) on local proc + std::complex **Rqsum; // Accumulator for conj(R(q)_alpha)*R(q)_beta + std::complex **Phi_q; // Phi's on local proc + std::complex **Phi_all; // Phi for all + + void readmap(); // to read the mapping of gf atoms + char *mapfile; // file name of the map file + + void getmass(); // to get the mass of each atom in a unit cell + + int nasr; + void postprocess(); // to post process the data + void EnforceASR(); // to apply acoustic sum rule to gamma point force constant matrix + + char *id_temp; // compute id for temperature + double *TempSum; // to get the average temperature vector + double inv_nTemp; // inverse of number of atoms in temperature group + class Compute *temperature; // compute that computes the temperature + + double hsum[6], **basis; + int *basetype; + + // private methods to do matrix inversion + void GaussJordan(int, std::complex*); + +}; +} +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal fix phonon command... + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: No atom found for fix phonon! + +Self-explanatory. Number of atoms in the group that was passed to +fix-phonon is less than 1. + +E: Can not open output file %s" + +Self-explanatory. + +E: Illegal fix_modify command + +Self-explanatory. + +E: Could not find fix_modify temp ID + +Self-explanatory. + +E: Fix_modify temp ID does not compute temperature + +Self-explanatory. + +E: Cannot open input map file %s + +Self-explanatory. + +E: Error while reading header of mapping file! + +Self-explanatory. The first line of the map file is expected to +contain 4 positive integer numbers. + +E: FFT mesh and number of atoms in group mismatch! + +Self-explanatory. The product of the 4 numbers should be exactly the +total number of atoms in the group that was passed to fix-phonon. + +E: Error while reading comment of mapping file! + +Self-explanatory. The second line of the map file should be a comment line. + +E: The mapping is incomplete! + +Self-explanatory. + +E: Error while reading mapping file! + +Self-explanatory. + +E: The mapping info read is incorrect! + +Self-explanatory. + +E: Singular matrix in complex GaussJordan! + +Self-explanatory. + +W: More than one fix phonon defined + +Self-explanatory. Just to warn that more than one fix-phonon is defined, but allowed. + +*/ diff --git a/src/dump_cfg.cpp b/src/dump_cfg.cpp index e5adf3e252..05f33997c9 100755 --- a/src/dump_cfg.cpp +++ b/src/dump_cfg.cpp @@ -13,6 +13,7 @@ /* ---------------------------------------------------------------------- Contributing author: Liang Wan (Chinese Academy of Sciences) + Memory efficiency improved by Ray Shan (Sandia) ------------------------------------------------------------------------- */ #include "math.h" @@ -34,7 +35,7 @@ using namespace LAMMPS_NS; -enum{INT,DOUBLE}; // same as in dump_custom.cpp +enum{INT,DOUBLE,STRING}; // same as in DumpCustom /* ---------------------------------------------------------------------- */ @@ -42,26 +43,23 @@ DumpCFG::DumpCFG(LAMMPS *lmp, int narg, char **arg) : DumpCustom(lmp, narg, arg) { if (narg < 10 || - strcmp(arg[5],"id") != 0 || strcmp(arg[6],"type") != 0 || + strcmp(arg[5],"mass") != 0 || strcmp(arg[6],"type") != 0 || (strcmp(arg[7],"xs") != 0 && strcmp(arg[7],"xsu") != 0) || (strcmp(arg[8],"ys") != 0 && strcmp(arg[8],"ysu") != 0) || (strcmp(arg[9],"zs") != 0 && strcmp(arg[9],"zsu") != 0)) error->all(FLERR,"Dump cfg arguments must start with " - "'id type xs ys zs' or 'id type xsu ysu zsu'"); + "'mass type xs ys zs' or 'mass type xsu ysu zsu'"); if (strcmp(arg[7],"xs") == 0) if (strcmp(arg[8],"ysu") == 0 || strcmp(arg[9],"zsu") == 0) - error->all(FLERR,"Dump cfg arguments can not mix xs|ys|zs with xsu|ysu|zsu"); + error->all(FLERR, + "Dump cfg arguments can not mix xs|ys|zs with xsu|ysu|zsu"); else unwrapflag = 0; else if (strcmp(arg[8],"ys") == 0 || strcmp(arg[9],"zs") == 0) - error->all(FLERR,"Dump cfg arguments can not mix xs|ys|zs with xsu|ysu|zsu"); + error->all(FLERR, + "Dump cfg arguments can not mix xs|ys|zs with xsu|ysu|zsu"); else unwrapflag = 1; - // arrays for data rearrangement - - rbuf = NULL; - nchosen = nlines = 0; - // setup auxiliary property name strings // convert 'X_ID[m]' (X=c,f,v) to 'ID_m' @@ -105,8 +103,6 @@ DumpCFG::DumpCFG(LAMMPS *lmp, int narg, char **arg) : DumpCFG::~DumpCFG() { - if (rbuf) memory->destroy(rbuf); - if (auxname) { for (int i = 0; i < nfield-5; i++) delete [] auxname[i]; delete [] auxname; @@ -117,7 +113,8 @@ DumpCFG::~DumpCFG() void DumpCFG::init_style() { - if (multifile == 0) error->all(FLERR,"Dump cfg requires one snapshot per file"); + if (multifile == 0) + error->all(FLERR,"Dump cfg requires one snapshot per file"); DumpCustom::init_style(); } @@ -160,78 +157,57 @@ void DumpCFG::write_header(bigint n) fprintf(fp,"entry_count = %d\n",nfield-2); for (int i = 0; i < nfield-5; i++) fprintf(fp,"auxiliary[%d] = %s\n",i,auxname[i]); - - // allocate memory needed for data rearrangement - - nchosen = static_cast (n); - if (rbuf) memory->destroy(rbuf); - memory->create(rbuf,nchosen,size_one,"dump:rbuf"); } -/* ---------------------------------------------------------------------- - write data lines to file in a block-by-block style - write head of block (mass & element name) only if has atoms of the type -------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- */ void DumpCFG::write_data(int n, double *mybuf) { - int i,j,m,itype; + int i,j,m; - double *rmass = atom->rmass; - double *mass = atom->mass; - - // transfer data from buf to rbuf - // if write by proc 0, transfer chunk by chunk - - for (i = 0, m = 0; i < n; i++) { - for (j = 0; j < size_one; j++) - rbuf[nlines][j] = mybuf[m++]; - nlines++; - } - - // write data lines in rbuf to file after transfer is done - - double unwrap_coord; - - if (nlines == nchosen) { - for (itype = 1; itype <= ntypes; itype++) { - for (i = 0; i < nchosen; i++) - if (rbuf[i][1] == itype) break; - if (i < nchosen) { - if (rmass) fprintf(fp,"%g\n",rmass[i]); - else fprintf(fp,"%g\n",mass[itype]); - fprintf(fp,"%s\n",typenames[itype]); - for (; i < nchosen; i++) { - if (rbuf[i][1] == itype) { - if (unwrapflag == 0) - for (j = 2; j < size_one; j++) { - if (vtype[j] == INT) - fprintf(fp,vformat[j],static_cast (rbuf[i][j])); - else fprintf(fp,vformat[j],rbuf[i][j]); - } - else { - - // Unwrapped scaled coordinates are shifted to - // center of expanded box, to prevent - // rewrapping by AtomEye. Dividing by - // expansion factor restores correct - // interatomic distances. - - for (j = 2; j < 5; j++) { - unwrap_coord = (rbuf[i][j] - 0.5)/UNWRAPEXPAND + 0.5; - fprintf(fp,vformat[j],unwrap_coord); - } - for (j = 5; j < size_one; j++) { - if (vtype[j] == INT) - fprintf(fp,vformat[j],static_cast (rbuf[i][j])); - else fprintf(fp,vformat[j],rbuf[i][j]); - } - } - fprintf(fp,"\n"); - } + if (unwrapflag == 0) { + m = 0; + for (i = 0; i < n; i++) { + for (j = 0; j < size_one; j++) { + if (j == 0) { + fprintf(fp,"%f \n",mybuf[m]); + } else if (j == 1) { + fprintf(fp,"%s \n",typenames[(int) mybuf[m]]); + } else if (j >= 2) { + if (vtype[j] == INT) + fprintf(fp,vformat[j],static_cast (mybuf[m])); + else if (vtype[j] == DOUBLE) + fprintf(fp,vformat[j],mybuf[m]); + else if (vtype[j] == STRING) + fprintf(fp,vformat[j],typenames[(int) mybuf[m]]); } + m++; } + fprintf(fp,"\n"); + } + } else if (unwrapflag == 1) { + m = 0; + double unwrap_coord; + for (i = 0; i < n; i++) { + for (j = 0; j < size_one; j++) { + if (j == 0) { + fprintf(fp,"%f \n",mybuf[m]); + } else if (j == 1) { + fprintf(fp,"%s \n",typenames[(int) mybuf[m]]); + } else if (j >= 2 && j <= 4) { + unwrap_coord = (mybuf[m] - 0.5)/UNWRAPEXPAND + 0.5; + fprintf(fp,vformat[j],unwrap_coord); + } else if (j >= 5 ) { + if (vtype[j] == INT) + fprintf(fp,vformat[j],static_cast (mybuf[m])); + else if (vtype[j] == DOUBLE) + fprintf(fp,vformat[j],mybuf[m]); + else if (vtype[j] == STRING) + fprintf(fp,vformat[j],typenames[(int) mybuf[m]]); + } + m++; + } + fprintf(fp,"\n"); } - nlines = 0; } } diff --git a/src/dump_cfg.h b/src/dump_cfg.h index 654582b1c6..d686b23287 100755 --- a/src/dump_cfg.h +++ b/src/dump_cfg.h @@ -31,14 +31,12 @@ class DumpCFG : public DumpCustom { private: char **auxname; // name strings of auxiliary properties - int nchosen; // # of lines to be written on a writing proc - int nlines; // # of lines transferred from buf to rbuf - double **rbuf; // buf of data lines for data lines rearrangement int unwrapflag; // 1 if unwrapped coordinates are requested void init_style(); void write_header(bigint); void write_data(int, double *); + }; } diff --git a/src/dump_custom.cpp b/src/dump_custom.cpp index 4e4b4224e4..e3682bbaca 100644 --- a/src/dump_custom.cpp +++ b/src/dump_custom.cpp @@ -44,7 +44,7 @@ enum{ID,MOL,TYPE,ELEMENT,MASS, TQX,TQY,TQZ,SPIN,ERADIUS,ERVEL,ERFORCE, COMPUTE,FIX,VARIABLE}; enum{LT,LE,GT,GE,EQ,NEQ}; -enum{INT,DOUBLE,STRING}; +enum{INT,DOUBLE,STRING}; // same as in DumpCFG #define INVOKED_PERATOM 8 diff --git a/tools/amber2lmp/README b/tools/amber2lmp/README index cf43a171a6..0c712f54ba 100644 --- a/tools/amber2lmp/README +++ b/tools/amber2lmp/README @@ -8,6 +8,14 @@ dump2try99.py same as dump2trj for LAMMPS 99 format ------------------------- +Modifications in file amber2lammps.py, by Vikas Varshney +Dated Nov 4, 2013 +Email address: vv0210@gmail.com + +added support for flags used in current version of AMBER + +------------------------- + Modifications in file amber2lammps.py, by Vikas Varshney Dated July 5, 2005 Email address : vv5@uakron.edu diff --git a/tools/amber2lmp/amber2lammps.py b/tools/amber2lmp/amber2lammps.py index 08bdf38d9e..cc74cf2bb2 100644 --- a/tools/amber2lmp/amber2lammps.py +++ b/tools/amber2lmp/amber2lammps.py @@ -8,6 +8,7 @@ # # Modified by Vikas Varshney, U Akron, 5 July 2005, as described in README # Bug Fixed :Third argument in Dihedral Coeffs section is an integer - Ketan S Khare September 26, 2011 +# Modified by Vikas Varshney, Oct 8, 2013 to include additional flags (Atomic_Number, Coulombic and van der Waals 1-4 factors which are included in newer vesions of .top and .crd files in amber12. #============================================================ @@ -554,6 +555,11 @@ class Amber: for i in range(self.NATOM): self.CHRG.append(eval(Pop(Item_list,0))) + print 'Reading Atomic Number...' + self.ANUMBER = [] + for i in range(self.NATOM): + self.ANUMBER.append(eval(Pop(Item_list,0))) + print 'Reading Atomic Masses...' self.AMASS = [] for i in range(self.NATOM): @@ -619,6 +625,16 @@ class Amber: for i in range(self.NPTRA): self.PHASE.append(eval(Pop(Item_list,0))) + print 'Reading 1-4 Electrostatic Scaling Factor...' + self.SCEEFAC = [] + for i in range(self.NPTRA): + self.SCEEFAC.append(eval(Pop(Item_list,0))) + + print 'Reading 1-4 Van der Waals Scaling Factor...' + self.SCNBFAC = [] + for i in range(self.NPTRA): + self.SCNBFAC.append(eval(Pop(Item_list,0))) + print 'Reading Solty...' #I think this is currently not used in AMBER. Check it out, though self.SOLTY = [] for i in range(self.NATYP): diff --git a/tools/msi2lmp/README b/tools/msi2lmp/README index ebb52f6bdd..a70fc3471b 100644 --- a/tools/msi2lmp/README +++ b/tools/msi2lmp/README @@ -1,6 +1,13 @@ Stephanie Teich-McGoldrick (Sandai) is the current maintainer of the msi2lmp tool. She can be contacted at steichm at sandia.gov +08 Oct 2013 Axel Kohlmeyer + +Fixed a memory access violation with Class 2 force fields. +Free all allocated memory to better detection of memory errors. +Print out version number and data with all print levels > 0. +Added valgrind checks to the regression tests + 02 Aug 2013 Axel Kohlmeyer Added rudimentary support for OPLS-AA based on diff --git a/tools/msi2lmp/TriclinicModification.pdf b/tools/msi2lmp/TriclinicModification.pdf deleted file mode 100644 index a88bb902f0..0000000000 Binary files a/tools/msi2lmp/TriclinicModification.pdf and /dev/null differ