From a0d74ca2ae7dcfdc75c1eec58c43607a4b203573 Mon Sep 17 00:00:00 2001 From: Michael Brown Date: Sun, 27 Oct 2019 22:31:00 -0700 Subject: [PATCH] Bug fix for gay-berne potential when mu != 1.0. --- doc/src/PDF/pair_gayberne_extra.pdf | Bin 99860 -> 99902 bytes doc/src/PDF/pair_gayberne_extra.tex | 2 +- lib/gpu/lal_gayberne.cu | 7 ++-- src/ASPHERE/pair_gayberne.cpp | 46 ++++++++++++------------- src/USER-INTEL/pair_gayberne_intel.cpp | 21 ++++++----- 5 files changed, 36 insertions(+), 40 deletions(-) diff --git a/doc/src/PDF/pair_gayberne_extra.pdf b/doc/src/PDF/pair_gayberne_extra.pdf index b794ac3382f2ed43d801ff97f68c81ae60eb8ac3..5ac8bd051bd371b5b6580ea8b67f03a0f3662afe 100644 GIT binary patch delta 2643 zcmV-Z3as^%j0V1p2C$?M0ys94pa>{`y;|Fj+eQ+7pReHCN{H3VzIVqkkTr4vjJzx| z0Tyt)fSi%z!DeP+jl7fn`#nWc6uVg@Tejkt(XFSex-O@xPV?pU)!)CDv650mYgKjq zu2L=v6;kYIX$o4~@8zCKM#H;VFuRdH!C{Zw_JS3$I-6aQZQboI}xm*cpGg2G7b<$?%*!3>q$bJkG` zY2d0LR4vo^R-^cKCsQ9IHx*7P74o3cNl>|ZA_o$~QM~uk*hu$6hgfaWtt59siu7J7 z!>MgWuzeV@(W`~B^0=br^ao6lM zD{lw{h+*L4z|hy6Qc*;I<0kx(Bae2y%zv+~lq~)@z=bUOV{u*l%<`c z?+|}`f5TQIL)C%S0(j9_WGR*I?1i!6jE3SRWBe#A@~Etza7;hpO8W_Yv3|mo_7k$C zpTORYy}5qk@UZ{zao1M#w$Qt_(EF~?yROiG&FSEr&`6Jsmgo(8S8Br^s}4Fri(5`z z7+fbaH|FxmMi3t!Z}G7f0r`p$s| zUrm_KJ2?jBOr+Km+2Q81V*uYaFYUTCl=Ms}Aw*Q*JkJJy6hhe5Zqf4^W0#Xw9` zKj(Q3JU2dmInTr8%Af7{$po||uqd07k+^|Bpo)uL_hYYT?!gqH;q?;tNa13fbx-h! ziSXszh0Ce&&hAY!zxG!~}J8JEq8eg+$pIZbiH@((`rq;*`a%sXbv_oQ|eOm+ht zR5K&a^x!&Z(v@NRD1YHOXrTQMyuFAq`XqZP)Q0TQi1>B?e6Tv#5FTsKIc^0 zNa=(@tFqBS=sJ|Bre{K(7gdBW$!^{jC93z!lMO1EPuY zjGw!dglRP~ZLaab+}vde#m?pAOaUk7&`sh4`7YKi|UqdAuN!9sS>e=<8bGq+AwEx$bB!9+PGUalTZ@pAsJsz)yw zfkZfbUO5`Iv63r$HXf|_Un#~Sr zwh1aHqBi*0zfp_P;Ff4DR~?(pmJ5FE1E+d_qFJrSh`cW%}JUYQ1ki>>udU)9`X{uyp@$*dF>Lcc$X* z;AEs~8#`m!6t{8{dK+(M4 z@W@^JeaC7~U63Ueqm5}6Wk{?8n?)mTy{=7m2EI^Yr8kJ^I)qXRZaLr@)CDk?1RN4x znpsxjxY=0&cb@q&5TXz$3^`lc;4OumD+{?$5^}yYD{zRVQY?h}n_+fJ*IEci5-*RA{-@LW>|nu3uWr>X7NjZA6zyv08#Rx+v?dvENbg z-nB*dak}RDrsjpV=J~$nZAVYMoEn+g+DFV^I6Y1c`po^JWM~SF)6|=ow&_8Gv*6J!4%{dHq&F z;S@8IXG7Zgl;$XS8f5Ag2LNDf0{!YyK5}U8&=0v#a8`phFS-oO(`SN?dkqeYUUN!E zpjeQyl*wGa-cs51R(oe{(Qbv*uIvoc&ezGjt1VK?_2HJzv{IT`DgmZ{r~4BAL8I)E zLJ_g+z*mjvG4*C@;kCPntNOjpaMhxx)TdrVO5T8`wtPa&P%+67S zFfceVGA%GSConK4DGD!5Z)8MabY&nYL^?7sF*z|XGB`0XFgP(Xx7B$8{t8kxMM60; zMm0G^LO3usLqIv%0VP%%qqW`0N}c|F4|i9qC3kkUiru_eH<~k?%b7F7eRh8K&5tTpN{MK#uFl_8 z+C`zGuAGh{P^iu?s#lBOq&Caf=g+a&^%pBCq^o4J=t;9e+iZmv2H&(xerlIcf8|B= zQnh!fTMHLNg^s+AHRFPX^PkT{oiP~x!eTh4?lEn~d|Z$#AXDbC+-mMhY3ak1H6{1hnby#e!n~!gn z`13C(Z{;8_Se$@_H?4iYTdtsWPqy#-+z913BcAoUEs&(l@t_1byqLFsvugW6aQ7cJp1|V*)Xo5pfC!1xgbI?eI>V?e|5A*8n`M5 zEma=hdJy04VCsG3=E74-g)*op3o0i^av(EYiuX<$8|hw9h}A0JN^vJlk=`j~IJK<^ zwhJRRdP!)j_A4Z(PY{O%Q<`)+_lnab$xG;yLFqRB9>D5E}MK)1@W4l}~C^7=r4 z7zRG}41LZi6-6{|#2+Q{e~{~B{=T$IvH0VV+p$Ru^|BFhMu}e+quh z83sYM_O9*C z@rm8d_QS_bQ_-73@0vpI+d}W!LjNtNgL8tA9xGDdhMgkz^E?|+2w^APQO^_Fe=euyd0va){FLVjcy4_B za-N6FwLjVMqX}r|z@n;4M&<_ofI1%ax*vMIa1Xi&g4gG`M+q0>q zM70eNFDsV7%0F-))YJYa-d<^9^jY@Otc$Xq@EzalEE_O61_-RHOk)1fvW6>Qee^3R`zSxd!D#r(nBRY{FzI^eV zLHM)~rlRQU`R4>Cp}KlMjqu>$N~N*^=ry5tgcbS4S2eH#uF@E;H)u+E#`j%X!L*bn zZLapgRNob%Fxc_IhjLmWAQcQmv8H+V8OmuCNdh%RNfj2#+^+)0dDIjkpJDb&`@qk? za8iAMf1=R=nEsprSUlP!e);*}1P-58w+1#UQE{{)(MhA!iOHP*zbeu*MqnB_+^--F zNNqNx4Pdc;ziwWx z*1PTHr}SPgHlP09Zf;jCv9$<+nooe!n>G|~FK^m$$i>HQbke>a?q-sTm)Gm-q|#0- zSHR}Qi_MO;fY+eSU3jk<&YVba-tt`WAO>#jvtIX20wWWKdy_uh_3w~+ ze@nMFyY;R+a%(E?sz=3y=`lCfQ|2(?FhHS6Y@m=GrM2Su`_1(umJy_iy1#2OeB+&$ zwbXdq8Fx+8ToGEV563z?>-DkAreIJjjjE@6U!w6g)bH$OjM0ukI?X~d8dSEVhQqF1 zT}z+aWxp~q5m2eEtjCot%e!K=Pqo%2e-)#Rc@`C13j7(a7cJ%CW+L3vy%euH1ic8ltN%;$Q|!~{k)KK^Fl7n2{}JELV3T7hus<}c+eI0Ni5 z0dNk#7sI_=Ubj`yIK@PJQ_qHcf8I^cc<5m4vH$>JYzF=MR^4)F?!bp!m~mDEnipLG z=J}&1`@MRH1=k$Y5hxavnmaH%ev6`uTa@N%)Pin>+^+16bklESO%++Bmg~YTZ=94e zll>}8pT_e1^DM=aYa(LTp08@r!wHX(!jpRtS9N=x;HsmZ(w}?ZrFa7xe`)!In7FC_ zfgLG{phGdXEX#;T8A+4>!$lp!V^mf<&+fAlnOg#}g%eE~_UoS4%S8GkUH7vi)qD_7COMvY)W21Wpzoobi$dhf9qs zA)8Uz?svBdlWU1HrbdwFGTzgAp22;uQT@>&p8*DXr9m?xBtpxhXm5=F0MItA$Ah{V zx4Ibt8DI!8I59LeGBh}s+j#*S2r)P@G&M3bIJf?J0sab7Lqsz~FgZCiMm0e*L`5(* zMK?w;G(<*4GB+?WMlduoJ|H|pL^DJ%IXN^&H9<2(MKCr+H%2ftL`Fq2H!v|qFf=kg zT?#KuWo~D5Xdp2*F(8-0KmjO!9m_is!2lG6;d9~@gv48jMO|rc)9LmllQr9|<%$3MC~)PeuyQhvf$V diff --git a/doc/src/PDF/pair_gayberne_extra.tex b/doc/src/PDF/pair_gayberne_extra.tex index a63ab2b1ce..0568b54dcf 100644 --- a/doc/src/PDF/pair_gayberne_extra.tex +++ b/doc/src/PDF/pair_gayberne_extra.tex @@ -131,7 +131,7 @@ and $$ \frac{ \partial \chi_{12} }{ \partial \mathbf{q}_i } = 4.0 \cdot 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 expression due to the determinant. Let $a_{mi}$ be the mth row of the diff --git a/lib/gpu/lal_gayberne.cu b/lib/gpu/lal_gayberne.cu index dc6e00ec82..cd1ee59fc6 100644 --- a/lib/gpu/lal_gayberne.cu +++ b/lib/gpu/lal_gayberne.cu @@ -316,10 +316,9 @@ __kernel void k_gayberne(const __global numtyp4 *restrict x_, numtyp tempv[3]; gpu_row_times3(iota,b1,tempv); gpu_cross3(tempv,iota,tchi); - temp1 = (numtyp)-4.0*ir*ir; - tchi[0] *= temp1; - tchi[1] *= temp1; - tchi[2] *= temp1; + tchi[0] *= temp2; + tchi[1] *= temp2; + tchi[2] *= temp2; } numtyp temp2 = factor_lj*eta*chi; diff --git a/src/ASPHERE/pair_gayberne.cpp b/src/ASPHERE/pair_gayberne.cpp index ea5a7e9aaa..7dd5270309 100644 --- a/src/ASPHERE/pair_gayberne.cpp +++ b/src/ASPHERE/pair_gayberne.cpp @@ -28,7 +28,6 @@ #include "citeme.h" #include "memory.h" #include "error.h" -#include "utils.h" using namespace LAMMPS_NS; @@ -461,20 +460,20 @@ void PairGayBerne::read_restart(FILE *fp) int i,j; int me = comm->me; 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); 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); } 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); if (setflag[i][j]) { if (me == 0) { - utils::sfread(FLERR,&epsilon[i][j],sizeof(double),1,fp,NULL,error); - utils::sfread(FLERR,&sigma[i][j],sizeof(double),1,fp,NULL,error); - utils::sfread(FLERR,&cut[i][j],sizeof(double),1,fp,NULL,error); + fread(&epsilon[i][j],sizeof(double),1,fp); + fread(&sigma[i][j],sizeof(double),1,fp); + fread(&cut[i][j],sizeof(double),1,fp); } MPI_Bcast(&epsilon[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; if (me == 0) { - utils::sfread(FLERR,&gamma,sizeof(double),1,fp,NULL,error); - utils::sfread(FLERR,&upsilon,sizeof(double),1,fp,NULL,error); - utils::sfread(FLERR,&mu,sizeof(double),1,fp,NULL,error); - utils::sfread(FLERR,&cut_global,sizeof(double),1,fp,NULL,error); - utils::sfread(FLERR,&offset_flag,sizeof(int),1,fp,NULL,error); - utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,NULL,error); + fread(&gamma,sizeof(double),1,fp); + fread(&upsilon,sizeof(double),1,fp); + fread(&mu,sizeof(double),1,fp); + fread(&cut_global,sizeof(double),1,fp); + fread(&offset_flag,sizeof(int),1,fp); + fread(&mix_flag,sizeof(int),1,fp); } MPI_Bcast(&gamma,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]); temp1 = -eta*u_r; - temp2 = eta*chi; - fforce[0] = temp1*dchi[0]-temp2*dUr[0]; - fforce[1] = temp1*dchi[1]-temp2*dUr[1]; - fforce[2] = temp1*dchi[2]-temp2*dUr[2]; + temp3 = eta*chi; + fforce[0] = temp1*dchi[0]-temp3*dUr[0]; + fforce[1] = temp1*dchi[1]-temp3*dUr[1]; + fforce[2] = temp1*dchi[2]-temp3*dUr[2]; // torque for particle 1 and 2 // 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::cross3(tempv,iota,dchi); - temp1 = -4.0/rsq; - dchi[0] *= temp1; - dchi[1] *= temp1; - dchi[2] *= temp1; + dchi[0] *= temp2; + dchi[1] *= temp2; + dchi[2] *= temp2; double dchi2[3]; if (newton_pair || j < nlocal) { MathExtra::vecmat(iota,b2,tempv); MathExtra::cross3(tempv,iota,dchi2); - dchi2[0] *= temp1; - dchi2[1] *= temp1; - dchi2[2] *= temp1; + dchi2[0] *= temp2; + dchi2[1] *= temp2; + dchi2[2] *= temp2; } // compute d_eta diff --git a/src/USER-INTEL/pair_gayberne_intel.cpp b/src/USER-INTEL/pair_gayberne_intel.cpp index 862dee2287..30e61f67b9 100644 --- a/src/USER-INTEL/pair_gayberne_intel.cpp +++ b/src/USER-INTEL/pair_gayberne_intel.cpp @@ -555,10 +555,10 @@ void PairGayBerneIntel::eval(const int offload, const int vflag, dchi_2 = temp2 * (iota_2 - temp1 * r12hat_2); temp1 = -eta * u_r; - temp2 = eta * chi; - fforce_0 = temp1 * dchi_0 - temp2 * dUr_0; - fforce_1 = temp1 * dchi_1 - temp2 * dUr_1; - fforce_2 = temp1 * dchi_2 - temp2 * dUr_2; + temp3 = eta * chi; + fforce_0 = temp1 * dchi_0 - temp3 * dUr_0; + fforce_1 = temp1 * dchi_1 - temp3 * dUr_1; + fforce_2 = temp1 * dchi_2 - temp3 * dUr_2; // torque for particle 1 and 2 // compute dUr @@ -579,18 +579,17 @@ void PairGayBerneIntel::eval(const int offload, const int vflag, ME_vecmat(iota, b1, tempv); ME_cross3(tempv, iota, dchi); - temp1 = (flt_t)-4.0 / rsq_form[jj]; - dchi_0 *= temp1; - dchi_1 *= temp1; - dchi_2 *= temp1; + dchi_0 *= temp2; + dchi_1 *= temp2; + dchi_2 *= temp2; flt_t dchi2_0, dchi2_1, dchi2_2; if (NEWTON_PAIR) { ME_vecmat(iota, b2, tempv); ME_cross3(tempv, iota, dchi2); - dchi2_0 *= temp1; - dchi2_1 *= temp1; - dchi2_2 *= temp1; + dchi2_0 *= temp2; + dchi2_1 *= temp2; + dchi2_2 *= temp2; } // compute d_eta