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