Commit 44bda245 by Steven J. Plimpton

### resolution of cutoff issues

parent 59d3af44
 ... @@ -10,8 +10,6 @@ pair_style atm command :h3 ... @@ -10,8 +10,6 @@ pair_style atm command :h3 [Syntax:] [Syntax:] SJP: add an arg pair_style atm cutoff cutoff_triple :pre pair_style atm cutoff cutoff_triple :pre cutoff = cutoff for each pair in 3-body interaction (distance units) cutoff = cutoff for each pair in 3-body interaction (distance units) ... @@ -19,11 +17,10 @@ cutoff_triple = additional cutoff applied to product of 3 pairwise distances (di ... @@ -19,11 +17,10 @@ cutoff_triple = additional cutoff applied to product of 3 pairwise distances (di [Examples:] [Examples:] SJP: is 1.5 a good value? pair_style atm 4.5 2.5 pair_style atm 2.5 1.5 pair_coeff * * * 0.072 :pre pair_coeff * * * 0.072 :pre pair_style hybrid/overlay lj/cut 6.5 atm 2.5 1.5 pair_style hybrid/overlay lj/cut 6.5 atm 4.5 2.5 pair_coeff * * lj/cut 1.0 1.0 pair_coeff * * lj/cut 1.0 1.0 pair_coeff 1 1 atm 1 0.064 pair_coeff 1 1 atm 1 0.064 pair_coeff 1 1 atm 2 0.080 pair_coeff 1 1 atm 2 0.080 ... @@ -55,9 +52,9 @@ command as in the example above. ... @@ -55,9 +52,9 @@ command as in the example above. The potential for a triplet of atom is calculated only if all 3 The potential for a triplet of atom is calculated only if all 3 distances r12, r23, r31 between the 3 atoms satisfy rIJ < cutoff. distances r12, r23, r31 between the 3 atoms satisfy rIJ < cutoff. SJP: In addition, the product of the 3 distances r12*r23*r31 < In addition, the product of the 3 distances r12*r23*r31 < cutoff_triple^3 is required, which limits the contribution of the cutoff_triple^3 is required, which excludes from calculation the potential to ??? triplets with small contribution to the interaction. The following coefficients must be defined for each pair of atoms The following coefficients must be defined for each pair of atoms types via the "pair_coeff"_pair_coeff.html command as in the examples types via the "pair_coeff"_pair_coeff.html command as in the examples ... ...
 ... @@ -16,7 +16,7 @@ region box block 0 \${xx} 0 \${yy} 0 \${zz} ... @@ -16,7 +16,7 @@ region box block 0 \${xx} 0 \${yy} 0 \${zz} create_box 1 box create_box 1 box create_atoms 1 box create_atoms 1 box pair_style hybrid/overlay lj/cut 4.5 atm 2.5 pair_style hybrid/overlay lj/cut 4.5 atm 4.5 2.5 pair_coeff * * lj/cut 1.0 1.0 pair_coeff * * lj/cut 1.0 1.0 pair_coeff * * atm * 0.072 pair_coeff * * atm * 0.072 ... ...
 LAMMPS (29 Jun 2018) LAMMPS (22 Aug 2018) # Axilrod-Teller-Muto potential example # Axilrod-Teller-Muto potential example variable x index 1 variable x index 1 ... @@ -26,9 +26,9 @@ Created orthogonal box = (0 0 0) to (18.3252 18.3252 18.3252) ... @@ -26,9 +26,9 @@ Created orthogonal box = (0 0 0) to (18.3252 18.3252 18.3252) 1 by 1 by 1 MPI processor grid 1 by 1 by 1 MPI processor grid create_atoms 1 box create_atoms 1 box Created 4000 atoms Created 4000 atoms Time spent = 0.00125098 secs Time spent = 0.00126314 secs pair_style hybrid/overlay lj/cut 4.5 atm 2.5 pair_style hybrid/overlay lj/cut 4.5 atm 4.5 2.5 pair_coeff * * lj/cut 1.0 1.0 pair_coeff * * lj/cut 1.0 1.0 pair_coeff * * atm * 0.072 pair_coeff * * atm * 0.072 ... @@ -60,26 +60,26 @@ Neighbor list info ... ... @@ -60,26 +60,26 @@ Neighbor list info ... bin: standard bin: standard Per MPI rank memory allocation (min/avg/max) = 11.47 | 11.47 | 11.47 Mbytes Per MPI rank memory allocation (min/avg/max) = 11.47 | 11.47 | 11.47 Mbytes Step Temp E_pair E_mol TotEng Press Step Temp E_pair E_mol TotEng Press 0 1.033 -4.733482 0 -3.1843694 -3.924644 0 1.033 -4.8899813 0 -3.3408686 -4.2298176 5 1.0335743 -4.7330528 0 -3.1830789 -3.9119065 5 1.0337853 -4.8928208 0 -3.3425304 -4.2233154 10 1.0349987 -4.7346788 0 -3.1825689 -3.8769962 10 1.0358056 -4.8953304 0 -3.3420104 -4.1897183 15 1.0362024 -4.7401425 0 -3.1862275 -3.8225106 15 1.0380938 -4.8990457 0 -3.3422942 -4.1310148 20 1.0352365 -4.7459627 0 -3.1934962 -3.7403572 20 1.0389566 -4.9014345 0 -3.3433892 -4.0406616 25 1.0295963 -4.7444397 0 -3.2004313 -3.6132651 25 1.0358313 -4.8989663 0 -3.3456079 -3.9093019 Loop time of 27.841 on 1 procs for 25 steps with 4000 atoms Loop time of 12.2062 on 1 procs for 25 steps with 4000 atoms Performance: 155.167 tau/day, 0.898 timesteps/s Performance: 353.920 tau/day, 2.048 timesteps/s 100.0% CPU use with 1 MPI tasks x no OpenMP threads 99.9% CPU use with 1 MPI tasks x no OpenMP threads MPI task timing breakdown: MPI task timing breakdown: Section | min time | avg time | max time |%varavg| %total Section | min time | avg time | max time |%varavg| %total --------------------------------------------------------------- --------------------------------------------------------------- Pair | 27.837 | 27.837 | 27.837 | 0.0 | 99.99 Pair | 12.202 | 12.202 | 12.202 | 0.0 | 99.96 Neigh | 0 | 0 | 0 | 0.0 | 0.00 Neigh | 0 | 0 | 0 | 0.0 | 0.00 Comm | 0.0015321 | 0.0015321 | 0.0015321 | 0.0 | 0.01 Comm | 0.0015621 | 0.0015621 | 0.0015621 | 0.0 | 0.01 Output | 0.00016594 | 0.00016594 | 0.00016594 | 0.0 | 0.00 Output | 0.00020814 | 0.00020814 | 0.00020814 | 0.0 | 0.00 Modify | 0.001785 | 0.001785 | 0.001785 | 0.0 | 0.01 Modify | 0.0019698 | 0.0019698 | 0.0019698 | 0.0 | 0.02 Other | | 0.0006731 | | | 0.00 Other | | 0.0007734 | | | 0.01 Nlocal: 4000 ave 4000 max 4000 min Nlocal: 4000 ave 4000 max 4000 min Histogram: 1 0 0 0 0 0 0 0 0 0 Histogram: 1 0 0 0 0 0 0 0 0 0 ... @@ -97,4 +97,4 @@ Dangerous builds = 0 ... @@ -97,4 +97,4 @@ Dangerous builds = 0 Please see the log.cite file for references relevant to this simulation Please see the log.cite file for references relevant to this simulation Total wall time: 0:00:29 Total wall time: 0:00:13
 LAMMPS (29 Jun 2018) LAMMPS (22 Aug 2018) # Axilrod-Teller-Muto potential example # Axilrod-Teller-Muto potential example variable x index 1 variable x index 1 ... @@ -26,9 +26,9 @@ Created orthogonal box = (0 0 0) to (18.3252 18.3252 18.3252) ... @@ -26,9 +26,9 @@ Created orthogonal box = (0 0 0) to (18.3252 18.3252 18.3252) 1 by 2 by 2 MPI processor grid 1 by 2 by 2 MPI processor grid create_atoms 1 box create_atoms 1 box Created 4000 atoms Created 4000 atoms Time spent = 0.000769138 secs Time spent = 0.000785112 secs pair_style hybrid/overlay lj/cut 4.5 atm 2.5 pair_style hybrid/overlay lj/cut 4.5 atm 4.5 2.5 pair_coeff * * lj/cut 1.0 1.0 pair_coeff * * lj/cut 1.0 1.0 pair_coeff * * atm * 0.072 pair_coeff * * atm * 0.072 ... @@ -60,26 +60,26 @@ Neighbor list info ... ... @@ -60,26 +60,26 @@ Neighbor list info ... bin: standard bin: standard Per MPI rank memory allocation (min/avg/max) = 5.532 | 5.532 | 5.532 Mbytes Per MPI rank memory allocation (min/avg/max) = 5.532 | 5.532 | 5.532 Mbytes Step Temp E_pair E_mol TotEng Press Step Temp E_pair E_mol TotEng Press 0 1.033 -4.733482 0 -3.1843694 -3.924644 0 1.033 -4.8921547 0 -3.343042 -4.2340557 5 1.0335743 -4.7330528 0 -3.1830789 -3.9119065 5 1.0337949 -4.8947881 0 -3.3444835 -4.2271456 10 1.0349987 -4.7346788 0 -3.1825689 -3.8769962 10 1.0358286 -4.8973178 0 -3.3439632 -4.1935779 15 1.0362024 -4.7401425 0 -3.1862275 -3.8225106 15 1.0381322 -4.9010593 0 -3.3442503 -4.134913 20 1.0352365 -4.7459627 0 -3.1934962 -3.7403572 20 1.0390107 -4.9034854 0 -3.3453589 -4.0446162 25 1.0295963 -4.7444397 0 -3.2004313 -3.6132651 25 1.0358988 -4.9010506 0 -3.3475908 -3.9133006 Loop time of 7.22029 on 4 procs for 25 steps with 4000 atoms Loop time of 3.20632 on 4 procs for 25 steps with 4000 atoms Performance: 598.314 tau/day, 3.462 timesteps/s Performance: 1347.340 tau/day, 7.797 timesteps/s 100.0% CPU use with 4 MPI tasks x no OpenMP threads 100.0% CPU use with 4 MPI tasks x no OpenMP threads MPI task timing breakdown: MPI task timing breakdown: Section | min time | avg time | max time |%varavg| %total Section | min time | avg time | max time |%varavg| %total --------------------------------------------------------------- --------------------------------------------------------------- Pair | 7.1346 | 7.1684 | 7.1863 | 0.8 | 99.28 Pair | 3.1207 | 3.1553 | 3.1859 | 1.5 | 98.41 Neigh | 0 | 0 | 0 | 0.0 | 0.00 Neigh | 0 | 0 | 0 | 0.0 | 0.00 Comm | 0.032945 | 0.0509 | 0.084664 | 9.1 | 0.70 Comm | 0.019466 | 0.05009 | 0.084602 | 12.0 | 1.56 Output | 0.00010133 | 0.00011051 | 0.00013804 | 0.0 | 0.00 Output | 7.1049e-05 | 8.2076e-05 | 0.00011325 | 0.0 | 0.00 Modify | 0.00059557 | 0.00060683 | 0.00061274 | 0.0 | 0.01 Modify | 0.00056338 | 0.00057292 | 0.00058413 | 0.0 | 0.02 Other | | 0.000318 | | | 0.00 Other | | 0.0003092 | | | 0.01 Nlocal: 1000 ave 1000 max 1000 min Nlocal: 1000 ave 1000 max 1000 min Histogram: 4 0 0 0 0 0 0 0 0 0 Histogram: 4 0 0 0 0 0 0 0 0 0 ... @@ -97,4 +97,4 @@ Dangerous builds = 0 ... @@ -97,4 +97,4 @@ Dangerous builds = 0 Please see the log.cite file for references relevant to this simulation Please see the log.cite file for references relevant to this simulation Total wall time: 0:00:07 Total wall time: 0:00:03
 ... @@ -15,9 +15,7 @@ ... @@ -15,9 +15,7 @@ Contributing author: Sergey Lishchuk Contributing author: Sergey Lishchuk ------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */ #include #include #include "pair_atm.h" #include "pair_atm.h" #include "atom.h" #include "atom.h" #include "citeme.h" #include "citeme.h" ... @@ -92,7 +90,6 @@ void PairATM::compute(int eflag, int vflag) ... @@ -92,7 +90,6 @@ void PairATM::compute(int eflag, int vflag) int *type = atom->type; int *type = atom->type; double cutoff_squared = cut_global*cut_global; double cutoff_squared = cut_global*cut_global; // SJP: new cutoff double triple = cut_triple*cut_triple*cut_triple; double triple = cut_triple*cut_triple*cut_triple; double cutoff_triple_sixth = triple*triple; double cutoff_triple_sixth = triple*triple; ... @@ -142,15 +139,14 @@ void PairATM::compute(int eflag, int vflag) ... @@ -142,15 +139,14 @@ void PairATM::compute(int eflag, int vflag) rik2 = rik[0]*rik[0] + rik[1]*rik[1] + rik[2]*rik[2]; rik2 = rik[0]*rik[0] + rik[1]*rik[1] + rik[2]*rik[2]; if (rik2 > cutoff_squared) continue; if (rik2 > cutoff_squared) continue; // SJP: add this line? double r6 = rij2*rjk2*rik2; if (r6 > cutoff_triple_sixth) continue; if (rij2*rjk2*rik2 > cutoff_triple_sixth) continue; nu_local = nu[type[i]][type[j]][type[k]]; nu_local = nu[type[i]][type[j]][type[k]]; if (nu_local == 0.0) continue; if (nu_local == 0.0) continue; interaction_ddd(nu_local, interaction_ddd(nu_local, rij2,rik2,rjk2,rij,rik,rjk,fj,fk,eflag,evdwl); r6,rij2,rik2,rjk2,rij,rik,rjk,fj,fk,eflag,evdwl); f[i][0] -= fj[0] + fk[0]; f[i][0] -= fj[0] + fk[0]; f[i][1] -= fj[1] + fk[1]; f[i][1] -= fj[1] + fk[1]; ... @@ -200,14 +196,10 @@ void PairATM::allocate() ... @@ -200,14 +196,10 @@ void PairATM::allocate() void PairATM::settings(int narg, char **arg) void PairATM::settings(int narg, char **arg) { { // SJP: change to 2 args, require triple <= global if (narg != 2) error->all(FLERR,"Illegal pair_style command"); if (narg != 2) error->all(FLERR,"Illegal pair_style command"); cut_global = force->numeric(FLERR,arg[0]); cut_global = force->numeric(FLERR,arg[0]); cut_triple = force->numeric(FLERR,arg[1]); cut_triple = force->numeric(FLERR,arg[1]); if (cut_triple > cut_global) error->all(FLERR,"Illegal pair_style command"); } } /* ---------------------------------------------------------------------- /* ---------------------------------------------------------------------- ... @@ -324,7 +316,6 @@ void PairATM::read_restart(FILE *fp) ... @@ -324,7 +316,6 @@ void PairATM::read_restart(FILE *fp) void PairATM::write_restart_settings(FILE *fp) void PairATM::write_restart_settings(FILE *fp) { { fwrite(&cut_global,sizeof(double),1,fp); fwrite(&cut_global,sizeof(double),1,fp); // SJP: 2nd arg fwrite(&cut_triple,sizeof(double),1,fp); fwrite(&cut_triple,sizeof(double),1,fp); } } ... @@ -334,7 +325,6 @@ void PairATM::write_restart_settings(FILE *fp) ... @@ -334,7 +325,6 @@ void PairATM::write_restart_settings(FILE *fp) void PairATM::read_restart_settings(FILE *fp) void PairATM::read_restart_settings(FILE *fp) { { // SJP: 2nd arg int me = comm->me; int me = comm->me; if (me == 0) { if (me == 0) { fread(&cut_global,sizeof(double),1,fp); fread(&cut_global,sizeof(double),1,fp); ... @@ -348,13 +338,12 @@ void PairATM::read_restart_settings(FILE *fp) ... @@ -348,13 +338,12 @@ void PairATM::read_restart_settings(FILE *fp) Axilrod-Teller-Muto (dipole-dipole-dipole) potential Axilrod-Teller-Muto (dipole-dipole-dipole) potential ------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */ void PairATM::interaction_ddd(double nu, void PairATM::interaction_ddd(double nu, double r6, double rij2, double rik2, double rjk2, double rij2, double rik2, double rjk2, double *rij, double *rik, double *rjk, double *rij, double *rik, double *rjk, double *fj, double *fk, int eflag, double &eng) double *fj, double *fk, int eflag, double &eng) { { double r6,r5inv,rri,rrj,rrk,rrr; double r5inv,rri,rrj,rrk,rrr; r6 = rij2*rik2*rjk2; r5inv = nu / (r6*r6*sqrt(r6)); r5inv = nu / (r6*r6*sqrt(r6)); rri = rik[0]*rij[0] + rik[1]*rij[1] + rik[2]*rij[2]; rri = rik[0]*rij[0] + rik[1]*rij[1] + rik[2]*rij[2]; rrj = rij[0]*rjk[0] + rij[1]*rjk[1] + rij[2]*rjk[2]; rrj = rij[0]*rjk[0] + rij[1]*rjk[1] + rij[2]*rjk[2]; ... ...
 ... @@ -39,12 +39,11 @@ class PairATM : public Pair { ... @@ -39,12 +39,11 @@ class PairATM : public Pair { void read_restart_settings(FILE *); void read_restart_settings(FILE *); protected: protected: // SJP: add 2nd cutoff double cut_global,cut_triple; double cut_global,cut_triple; double ***nu; double ***nu; void allocate(); void allocate(); void interaction_ddd(double, double, double, double, double *, void interaction_ddd(double, double, double, double, double, double *, double *, double *, double *, double *, int, double &); double *, double *, double *, double *, int, double &); }; }; ... ...
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment