ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Added hbnewflag argument to pair_style reax

Browse Source 
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@4849 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
athomps
2010-09-24 22:54:26 +00:00
parent 38867b49df
commit 68a5dad07d
2 changed files with 23 additions and 9 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
lib/reax/reax_connect.F
View File

@ -186,7 +186,7 @@ c$$$* end if
**********************************************************************
**********************************************************************
subroutine srtbon1(lprune,lhb,hbcut_in)
subroutine srtbon1(lprune,lhb,hbcut_in,lhbnew_in)
**********************************************************************
#include "cbka.blk"
@ -230,8 +230,9 @@ c$$$ call timer(65)
c$$$ close (65)
c$$$ end if
c Transfer hbcut from C++ calling function
c Transfer hbcut and lhbnew from C++ calling function
hbcut = hbcut_in
lhbnew = lhbnew_in
do i1=1,na
abo(i1)=0.0d0
@ -1528,6 +1529,7 @@ c
ihb(nhb,8)=k3
* write (64,*)nhb,ihb(nhb,1),j3,j1,j2,nbohb,k1,k2,k3,bo(nbohb),
* $dishb
end if
10 continue

26
lib/reax/reax_poten.F
View File

@ -2473,7 +2473,11 @@ c$$$ end if
sin2=sinhu*sinhu
exphu1=exp(-vhb1(ityhb)*boa)
exphu2=exp(-vhb2(ityhb)*(rhu1+rhu2-2.0))
ehbh=(1.0-exphu1)*dehb(ityhb)*exphu2*sin2*sin2*sin2*sin2
if (lhbnew .eq. 0) then
ehbh=(1.0-exphu1)*dehb(ityhb)*exphu2*sin2*sin2*sin2*sin2
else
ehbh=(1.0-exphu1)*dehb(ityhb)*exphu2*sin2*sin2
endif
ehb=ehb+ehbh
estrain(j(2))=estrain(j(2))+ehbh !2nd atom energy
@ -2482,12 +2486,20 @@ c$$$ end if
* Calculate first derivatives *
* *
**********************************************************************
dehbdbo=vhb1(ityhb)*exphu1*dehb(ityhb)*exphu2*sin2*sin2*
$sin2*sin2
dehbdv=(1.0-exphu1)*dehb(ityhb)*exphu2*
$4.0*sin2*sin2*sin2*sinhu*cos(hhb(i1)/2.0)
dehbdrda=(1.0-exphu1)*dehb(ityhb)*sin2*sin2*sin2*sin2*
$vhb2(ityhb)*(rhb(ityhb)/(rda*rda)-1.0/rhb(ityhb))*exphu2
if (lhbnew .eq. 0) then
dehbdbo=vhb1(ityhb)*exphu1*dehb(ityhb)*exphu2*sin2*sin2*
$ sin2*sin2
dehbdv=(1.0-exphu1)*dehb(ityhb)*exphu2*
$ 4.0*sin2*sin2*sin2*sinhu*cos(hhb(i1)/2.0)
dehbdrda=(1.0-exphu1)*dehb(ityhb)*sin2*sin2*sin2*sin2*
$ vhb2(ityhb)*(rhb(ityhb)/(rda*rda)-1.0/rhb(ityhb))*exphu2
else
dehbdbo=vhb1(ityhb)*exphu1*dehb(ityhb)*exphu2*sin2*sin2
dehbdv=(1.0-exphu1)*dehb(ityhb)*exphu2*
$ 2.0*sin2*sinhu*cos(hhb(i1)/2.0)
dehbdrda=(1.0-exphu1)*dehb(ityhb)*sin2*sin2*
$ vhb2(ityhb)*(rhb(ityhb)/(rda*rda)-1.0/rhb(ityhb))*exphu2
endif
if (Lvirial.eq.1) then
do k1=1,3