jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Bug fix for gay-berne potential when mu != 1.0.

This commit is contained in:
Michael Brown 2019-10-27 22:31:00 -07:00
parent 67b174701e
commit a0d74ca2ae
5 changed files with 36 additions and 40 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
doc/src/PDF/pair_gayberne_extra.pdf
View File

Binary file not shown.

2
doc/src/PDF/pair_gayberne_extra.tex
View File

@ -131,7 +131,7 @@ and
$$ \frac{ \partial \chi_{12} }{ \partial \mathbf{q}_i } = 4.0 \cdot $$ \frac{ \partial \chi_{12} }{ \partial \mathbf{q}_i } = 4.0 \cdot
r^{-2} \cdot \mathbf{A}_i (- \mathbf{\iota}^T \cdot \mathbf{B}_i r^{-2} \cdot \mathbf{A}_i (- \mathbf{\iota}^T \cdot \mathbf{B}_i
\times \mathbf{\iota} ). $$ \times \mathbf{\iota} ) \cdot \mu \cdot \chi_{12}^{ ( \mu -1 ) / \mu}. $$
For the derivative of the $\eta$ term, we were unable to find a matrix For the derivative of the $\eta$ term, we were unable to find a matrix
expression due to the determinant. Let $a_{mi}$ be the mth row of the expression due to the determinant. Let $a_{mi}$ be the mth row of the

7
lib/gpu/lal_gayberne.cu
View File

@ -316,10 +316,9 @@ __kernel void k_gayberne(const __global numtyp4 *restrict x_,
numtyp tempv[3]; numtyp tempv[3];
gpu_row_times3(iota,b1,tempv); gpu_row_times3(iota,b1,tempv);
gpu_cross3(tempv,iota,tchi); gpu_cross3(tempv,iota,tchi);
temp1 = (numtyp)-4.0*ir*ir; tchi[0] *= temp2;
tchi[0] *= temp1; tchi[1] *= temp2;
tchi[1] *= temp1; tchi[2] *= temp2;
tchi[2] *= temp1;
} }
numtyp temp2 = factor_lj*eta*chi; numtyp temp2 = factor_lj*eta*chi;

46
src/ASPHERE/pair_gayberne.cpp
View File

@ -28,7 +28,6 @@
#include "citeme.h" #include "citeme.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
#include "utils.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
@ -461,20 +460,20 @@ void PairGayBerne::read_restart(FILE *fp)
int i,j; int i,j;
int me = comm->me; int me = comm->me;
for (i = 1; i <= atom->ntypes; i++) { for (i = 1; i <= atom->ntypes; i++) {
if (me == 0) utils::sfread(FLERR,&setwell[i],sizeof(int),1,fp,NULL,error); if (me == 0) fread(&setwell[i],sizeof(int),1,fp);
MPI_Bcast(&setwell[i],1,MPI_INT,0,world); MPI_Bcast(&setwell[i],1,MPI_INT,0,world);
if (setwell[i]) { if (setwell[i]) {
if (me == 0) utils::sfread(FLERR,&well[i][0],sizeof(double),3,fp,NULL,error); if (me == 0) fread(&well[i][0],sizeof(double),3,fp);
MPI_Bcast(&well[i][0],3,MPI_DOUBLE,0,world); MPI_Bcast(&well[i][0],3,MPI_DOUBLE,0,world);
} }
for (j = i; j <= atom->ntypes; j++) { for (j = i; j <= atom->ntypes; j++) {
if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,NULL,error); if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
if (setflag[i][j]) { if (setflag[i][j]) {
if (me == 0) { if (me == 0) {
utils::sfread(FLERR,&epsilon[i][j],sizeof(double),1,fp,NULL,error); fread(&epsilon[i][j],sizeof(double),1,fp);
utils::sfread(FLERR,&sigma[i][j],sizeof(double),1,fp,NULL,error); fread(&sigma[i][j],sizeof(double),1,fp);
utils::sfread(FLERR,&cut[i][j],sizeof(double),1,fp,NULL,error); fread(&cut[i][j],sizeof(double),1,fp);
} }
MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
@ -506,12 +505,12 @@ void PairGayBerne::read_restart_settings(FILE *fp)
{ {
int me = comm->me; int me = comm->me;
if (me == 0) { if (me == 0) {
utils::sfread(FLERR,&gamma,sizeof(double),1,fp,NULL,error); fread(&gamma,sizeof(double),1,fp);
utils::sfread(FLERR,&upsilon,sizeof(double),1,fp,NULL,error); fread(&upsilon,sizeof(double),1,fp);
utils::sfread(FLERR,&mu,sizeof(double),1,fp,NULL,error); fread(&mu,sizeof(double),1,fp);
utils::sfread(FLERR,&cut_global,sizeof(double),1,fp,NULL,error); fread(&cut_global,sizeof(double),1,fp);
utils::sfread(FLERR,&offset_flag,sizeof(int),1,fp,NULL,error); fread(&offset_flag,sizeof(int),1,fp);
utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,NULL,error); fread(&mix_flag,sizeof(int),1,fp);
} }
MPI_Bcast(&gamma,1,MPI_DOUBLE,0,world); MPI_Bcast(&gamma,1,MPI_DOUBLE,0,world);
MPI_Bcast(&upsilon,1,MPI_DOUBLE,0,world); MPI_Bcast(&upsilon,1,MPI_DOUBLE,0,world);
@ -644,10 +643,10 @@ double PairGayBerne::gayberne_analytic(const int i,const int j,double a1[3][3],
dchi[2] = temp2*(iota[2]-temp1*r12hat[2]); dchi[2] = temp2*(iota[2]-temp1*r12hat[2]);
temp1 = -eta*u_r; temp1 = -eta*u_r;
temp2 = eta*chi; temp3 = eta*chi;
fforce[0] = temp1*dchi[0]-temp2*dUr[0]; fforce[0] = temp1*dchi[0]-temp3*dUr[0];
fforce[1] = temp1*dchi[1]-temp2*dUr[1]; fforce[1] = temp1*dchi[1]-temp3*dUr[1];
fforce[2] = temp1*dchi[2]-temp2*dUr[2]; fforce[2] = temp1*dchi[2]-temp3*dUr[2];
// torque for particle 1 and 2 // torque for particle 1 and 2
// compute dUr // compute dUr
@ -668,18 +667,17 @@ double PairGayBerne::gayberne_analytic(const int i,const int j,double a1[3][3],
MathExtra::vecmat(iota,b1,tempv); MathExtra::vecmat(iota,b1,tempv);
MathExtra::cross3(tempv,iota,dchi); MathExtra::cross3(tempv,iota,dchi);
temp1 = -4.0/rsq; dchi[0] *= temp2;
dchi[0] *= temp1; dchi[1] *= temp2;
dchi[1] *= temp1; dchi[2] *= temp2;
dchi[2] *= temp1;
double dchi2[3]; double dchi2[3];
if (newton_pair || j < nlocal) { if (newton_pair || j < nlocal) {
MathExtra::vecmat(iota,b2,tempv); MathExtra::vecmat(iota,b2,tempv);
MathExtra::cross3(tempv,iota,dchi2); MathExtra::cross3(tempv,iota,dchi2);
dchi2[0] *= temp1; dchi2[0] *= temp2;
dchi2[1] *= temp1; dchi2[1] *= temp2;
dchi2[2] *= temp1; dchi2[2] *= temp2;
} }
// compute d_eta // compute d_eta

21
src/USER-INTEL/pair_gayberne_intel.cpp
View File

@ -555,10 +555,10 @@ void PairGayBerneIntel::eval(const int offload, const int vflag,
dchi_2 = temp2 * (iota_2 - temp1 * r12hat_2); dchi_2 = temp2 * (iota_2 - temp1 * r12hat_2);
temp1 = -eta * u_r; temp1 = -eta * u_r;
temp2 = eta * chi; temp3 = eta * chi;
fforce_0 = temp1 * dchi_0 - temp2 * dUr_0; fforce_0 = temp1 * dchi_0 - temp3 * dUr_0;
fforce_1 = temp1 * dchi_1 - temp2 * dUr_1; fforce_1 = temp1 * dchi_1 - temp3 * dUr_1;
fforce_2 = temp1 * dchi_2 - temp2 * dUr_2; fforce_2 = temp1 * dchi_2 - temp3 * dUr_2;
// torque for particle 1 and 2 // torque for particle 1 and 2
// compute dUr // compute dUr
@ -579,18 +579,17 @@ void PairGayBerneIntel::eval(const int offload, const int vflag,
ME_vecmat(iota, b1, tempv); ME_vecmat(iota, b1, tempv);
ME_cross3(tempv, iota, dchi); ME_cross3(tempv, iota, dchi);
temp1 = (flt_t)-4.0 / rsq_form[jj]; dchi_0 *= temp2;
dchi_0 *= temp1; dchi_1 *= temp2;
dchi_1 *= temp1; dchi_2 *= temp2;
dchi_2 *= temp1;
flt_t dchi2_0, dchi2_1, dchi2_2; flt_t dchi2_0, dchi2_1, dchi2_2;
if (NEWTON_PAIR) { if (NEWTON_PAIR) {
ME_vecmat(iota, b2, tempv); ME_vecmat(iota, b2, tempv);
ME_cross3(tempv, iota, dchi2); ME_cross3(tempv, iota, dchi2);
dchi2_0 *= temp1; dchi2_0 *= temp2;
dchi2_1 *= temp1; dchi2_1 *= temp2;
dchi2_2 *= temp1; dchi2_2 *= temp2;
} }
// compute d_eta // compute d_eta