1 // clang-format off
2 /* ----------------------------------------------------------------------
3    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
4    https://www.lammps.org/, Sandia National Laboratories
5    Steve Plimpton, sjplimp@sandia.gov
6 
7    Copyright (2003) Sandia Corporation.  Under the terms of Contract
8    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
9    certain rights in this software.  This software is distributed under
10    the GNU General Public License.
11 
12    See the README file in the top-level LAMMPS directory.
13 ------------------------------------------------------------------------- */
14 
15 #include "bond_nonlinear.h"
16 
17 #include <cmath>
18 #include <cstring>
19 #include "atom.h"
20 #include "neighbor.h"
21 #include "comm.h"
22 #include "force.h"
23 #include "memory.h"
24 #include "error.h"
25 
26 
27 using namespace LAMMPS_NS;
28 
29 /* ---------------------------------------------------------------------- */
30 
BondNonlinear(LAMMPS * lmp)31 BondNonlinear::BondNonlinear(LAMMPS *lmp) : Bond(lmp) {}
32 
33 /* ---------------------------------------------------------------------- */
34 
~BondNonlinear()35 BondNonlinear::~BondNonlinear()
36 {
37   if (allocated) {
38     memory->destroy(setflag);
39     memory->destroy(epsilon);
40     memory->destroy(r0);
41     memory->destroy(lamda);
42   }
43 }
44 
45 /* ---------------------------------------------------------------------- */
46 
compute(int eflag,int vflag)47 void BondNonlinear::compute(int eflag, int vflag)
48 {
49   int i1,i2,n,type;
50   double delx,dely,delz,ebond,fbond;
51   double rsq,r,dr,drsq,lamdasq,denom,denomsq;
52 
53   ebond = 0.0;
54   ev_init(eflag,vflag);
55 
56   double **x = atom->x;
57   double **f = atom->f;
58   int **bondlist = neighbor->bondlist;
59   int nbondlist = neighbor->nbondlist;
60   int nlocal = atom->nlocal;
61   int newton_bond = force->newton_bond;
62 
63   for (n = 0; n < nbondlist; n++) {
64     i1 = bondlist[n][0];
65     i2 = bondlist[n][1];
66     type = bondlist[n][2];
67 
68     delx = x[i1][0] - x[i2][0];
69     dely = x[i1][1] - x[i2][1];
70     delz = x[i1][2] - x[i2][2];
71 
72     rsq = delx*delx + dely*dely + delz*delz;
73     r = sqrt(rsq);
74     dr = r - r0[type];
75     drsq = dr*dr;
76     lamdasq = lamda[type]*lamda[type];
77     denom = lamdasq - drsq;
78     denomsq = denom*denom;
79 
80     // force & energy
81 
82     fbond = -epsilon[type]/r * 2.0*dr*lamdasq/denomsq;
83     if (eflag) ebond = epsilon[type] * drsq / denom;
84 
85     // apply force to each of 2 atoms
86 
87     if (newton_bond || i1 < nlocal) {
88       f[i1][0] += delx*fbond;
89       f[i1][1] += dely*fbond;
90       f[i1][2] += delz*fbond;
91     }
92 
93     if (newton_bond || i2 < nlocal) {
94       f[i2][0] -= delx*fbond;
95       f[i2][1] -= dely*fbond;
96       f[i2][2] -= delz*fbond;
97     }
98 
99     if (evflag) ev_tally(i1,i2,nlocal,newton_bond,ebond,fbond,delx,dely,delz);
100   }
101 }
102 
103 /* ---------------------------------------------------------------------- */
104 
allocate()105 void BondNonlinear::allocate()
106 {
107   allocated = 1;
108   int n = atom->nbondtypes;
109 
110   memory->create(epsilon,n+1,"bond:epsilon");
111   memory->create(r0,n+1,"bond:r0");
112   memory->create(lamda,n+1,"bond:lamda");
113   memory->create(setflag,n+1,"bond:setflag");
114   for (int i = 1; i <= n; i++) setflag[i] = 0;
115 }
116 
117 /* ----------------------------------------------------------------------
118    set coeffs for one type
119 ------------------------------------------------------------------------- */
120 
coeff(int narg,char ** arg)121 void BondNonlinear::coeff(int narg, char **arg)
122 {
123   if (narg != 4) error->all(FLERR,"Incorrect args for bond coefficients");
124   if (!allocated) allocate();
125 
126   int ilo,ihi;
127   utils::bounds(FLERR,arg[0],1,atom->nbondtypes,ilo,ihi,error);
128 
129   double epsilon_one = utils::numeric(FLERR,arg[1],false,lmp);
130   double r0_one = utils::numeric(FLERR,arg[2],false,lmp);
131   double lamda_one = utils::numeric(FLERR,arg[3],false,lmp);
132 
133   int count = 0;
134   for (int i = ilo; i <= ihi; i++) {
135     epsilon[i] = epsilon_one;
136     r0[i] = r0_one;
137     lamda[i] = lamda_one;
138     setflag[i] = 1;
139     count++;
140   }
141 
142   if (count == 0) error->all(FLERR,"Incorrect args for bond coefficients");
143 }
144 
145 /* ---------------------------------------------------------------------- */
146 
equilibrium_distance(int i)147 double BondNonlinear::equilibrium_distance(int i)
148 {
149   return r0[i];
150 }
151 
152 /* ----------------------------------------------------------------------
153    proc 0 writes to restart file
154 ------------------------------------------------------------------------- */
155 
write_restart(FILE * fp)156 void BondNonlinear::write_restart(FILE *fp)
157 {
158   fwrite(&epsilon[1],sizeof(double),atom->nbondtypes,fp);
159   fwrite(&r0[1],sizeof(double),atom->nbondtypes,fp);
160   fwrite(&lamda[1],sizeof(double),atom->nbondtypes,fp);
161 }
162 
163 /* ----------------------------------------------------------------------
164    proc 0 reads from restart file, bcasts
165 ------------------------------------------------------------------------- */
166 
read_restart(FILE * fp)167 void BondNonlinear::read_restart(FILE *fp)
168 {
169   allocate();
170 
171   if (comm->me == 0) {
172     utils::sfread(FLERR,&epsilon[1],sizeof(double),atom->nbondtypes,fp,nullptr,error);
173     utils::sfread(FLERR,&r0[1],sizeof(double),atom->nbondtypes,fp,nullptr,error);
174     utils::sfread(FLERR,&lamda[1],sizeof(double),atom->nbondtypes,fp,nullptr,error);
175   }
176   MPI_Bcast(&epsilon[1],atom->nbondtypes,MPI_DOUBLE,0,world);
177   MPI_Bcast(&r0[1],atom->nbondtypes,MPI_DOUBLE,0,world);
178   MPI_Bcast(&lamda[1],atom->nbondtypes,MPI_DOUBLE,0,world);
179 
180   for (int i = 1; i <= atom->nbondtypes; i++) setflag[i] = 1;
181 }
182 
183 /* ----------------------------------------------------------------------
184    proc 0 writes to data file
185 ------------------------------------------------------------------------- */
186 
write_data(FILE * fp)187 void BondNonlinear::write_data(FILE *fp)
188 {
189   for (int i = 1; i <= atom->nbondtypes; i++)
190     fprintf(fp,"%d %g %g %g\n",i,epsilon[i],r0[i],lamda[i]);
191 }
192 
193 /* ---------------------------------------------------------------------- */
194 
single(int type,double rsq,int,int,double & fforce)195 double BondNonlinear::single(int type, double rsq, int /*i*/, int /*j*/,
196                              double &fforce)
197 {
198   double r = sqrt(rsq);
199   double dr = r - r0[type];
200   double drsq = dr*dr;
201   double lamdasq = lamda[type]*lamda[type];
202   double denom = lamdasq - drsq;
203   double denomsq = denom*denom;
204   fforce = -epsilon[type]/r * 2.0*dr*lamdasq/denomsq;
205   return epsilon[type] * drsq / denom;
206 }
207 
208 /* ---------------------------------------------------------------------- */
209 
extract(const char * str,int & dim)210 void *BondNonlinear::extract(const char *str, int &dim)
211 {
212   dim = 1;
213   if (strcmp(str,"epsilon")==0) return (void*) epsilon;
214   if (strcmp(str,"r0")==0) return (void*) r0;
215   return nullptr;
216 }
217