forked from lijiext/lammps
3986 lines
122 KiB
Fortran
3986 lines
122 KiB
Fortran
**********************************************************************
|
|
* *
|
|
* REAXFF Reactive force field program *
|
|
* *
|
|
* Developed and written by Adri van Duin, duin@wag.caltech.edu *
|
|
* *
|
|
* Copyright (c) 2001-2010 California Institute of Technology *
|
|
* *
|
|
* This is an open-source program. Feel free to modify its *
|
|
* contents. Please keep me informed of any useful modification *
|
|
* or addition that you made. Please do not distribute this *
|
|
* program to others; if people are interested in obtaining *
|
|
* a copy of this program let them contact me first. *
|
|
* *
|
|
**********************************************************************
|
|
********************************************************************
|
|
|
|
subroutine calval
|
|
|
|
**********************************************************************
|
|
#include "cbka.blk"
|
|
#include "cbkc.blk"
|
|
#include "cbkdhdc.blk"
|
|
#include "cbkdrdc.blk"
|
|
#include "cbkh.blk"
|
|
#include "cbkrbo.blk"
|
|
#include "cbkvalence.blk"
|
|
#include "cellcoord.blk"
|
|
#include "control.blk"
|
|
dimension a(3),b(3),j(3),dradc(3,3),drbdc(3,3),dtdc(3,3),
|
|
$dargdc(3,3),dndc(3,3),dadc(3),dbdc(3)
|
|
**********************************************************************
|
|
* *
|
|
* Calculate valency angles and their derivatives to cartesian *
|
|
* coordinates *
|
|
* Valency angle energies are calculated in valang *
|
|
* *
|
|
**********************************************************************
|
|
**********************************************************************
|
|
* Description of variables used in this routine.
|
|
*
|
|
* ndebug: stored in cbka.blk; control-parameter
|
|
* third: local variable
|
|
* twothird: local variable
|
|
* dadc(3): local array; stores derivative distance to cartesians
|
|
* dbdc(3): local array; stores derivative distance to cartesians
|
|
* i1: local do-loop counter
|
|
* i2: local do-loop counter
|
|
* k1: local do-loop counter
|
|
* k2: local do-loop counter
|
|
* dradc(3,3): local array; stores derivatives bond lengths to
|
|
* cartesians
|
|
* drbdc(3,3): local array; stores derivatives bond lengths to
|
|
* cartesians
|
|
* nval: stored in cbka.blk; number of valence angles
|
|
* ity: local integer; atom type
|
|
* iv(nvalmax,6): stored in cbka.blk; valence angle identifiers
|
|
* j(3): local integer array; stores valence angle atom numbers
|
|
* la: local integer: stores bond numbers in valence angle
|
|
* lb: local integer: stores bond numbers in valence angle
|
|
* ivl1: local integer; stores symmetric copy number of bond
|
|
* ivl2: local integer; stores symmetric copy number of bond
|
|
* ibsym(nbomax): stored in cbka.blk; symmetric copy number of bond
|
|
* isign1: local integer; -1 or 1
|
|
* isign2: local integer; -1 or 1
|
|
* rla: local variable; stores bond length for bond la
|
|
* rlb: local variable; stores bond length for bond lb
|
|
* rbo(nbomax): stored in cbka.blk; stores bond lengths
|
|
* ix1,iy1,iz1,ix2,iy2,iz2: local integers; periodic cell shifts
|
|
* a(3): local variable; distance in x,y and z-direction between atoms
|
|
* b(3): local variable; distance in x,y and z-direction between atoms
|
|
* c(nat,3): stored in cbka.blk; cartesian coordinate array
|
|
* tm11,tm21,tm22,tm31,tm32,tm33: stored in cbka.blk; periodic cell
|
|
* matrix
|
|
* poem: local variable; product of bond lengths
|
|
* tel: local variable; cross-product of x,y and z-interatomic
|
|
* distances
|
|
* arg: local variable; cosine of angle between bonds a and b
|
|
* arg2: local variable; square of arg
|
|
* s1ma22: local variable; used to check whether angle gets to 180
|
|
* degrees
|
|
* s1ma2: local variable; square root of s1ma22
|
|
* hl: local variable; angle (in radians) between bonds a and b
|
|
* h(nvamax): stored in cbka.blk; angle (in radians) between bonds a
|
|
* and b
|
|
* ib(nbomax,3): stored in cbka.blk: bond distance identifiers
|
|
* drdc(nbomax,3,2): stored in cbka.blk; derivatives bond distances
|
|
* to cartesian coordinates
|
|
* dndc(3,3): local variable; temporary storage for calculating
|
|
* derivatives of valence angle to cartesians
|
|
* dtdc(3,3): local variable; temporary storage for calculating
|
|
* derivatives of valence angle to cartesians
|
|
* dargdc(3,3): local variable; temporary storage for calculating
|
|
* derivatives of valence angle to cartesians
|
|
* dhdc(nvamax,3,3): stored in cbka.blk; derivatives of valence angle
|
|
* to cartesians
|
|
*
|
|
**********************************************************************
|
|
c$$$ if (ndebug.eq.1) then
|
|
c$$$C open (65,file='fort.65',status='unknown',access='append')
|
|
c$$$ write (65,*) 'In calval'
|
|
c$$$ call timer(65)
|
|
c$$$ close (65)
|
|
c$$$ end if
|
|
|
|
third=1.0/3.0
|
|
twothird=2.0/3.0
|
|
dadc(1)=-1.0
|
|
dadc(2)=1.0
|
|
dadc(3)=0.0
|
|
dbdc(1)=0.0
|
|
dbdc(2)=1.0
|
|
dbdc(3)=-1.0
|
|
do k1=1,3
|
|
do k2=1,3
|
|
dradc(k1,k2)=0.0
|
|
drbdc(k1,k2)=0.0
|
|
end do
|
|
end do
|
|
if (nval.eq.0) return
|
|
|
|
do 10 i1=1,nval
|
|
ity=iv(i1,1)
|
|
j(1)=iv(i1,2)
|
|
j(2)=iv(i1,3)
|
|
j(3)=iv(i1,4)
|
|
**********************************************************************
|
|
* *
|
|
* Determine valency angle *
|
|
* *
|
|
**********************************************************************
|
|
la=iv(i1,5)
|
|
lb=iv(i1,6)
|
|
ivl1=ibsym(la)
|
|
ivl2=ibsym(lb)
|
|
isign1=1
|
|
isign2=1
|
|
rla=rbo(la)
|
|
rlb=rbo(lb)
|
|
|
|
call dista2(j(2),j(1),dis,a(1),a(2),a(3))
|
|
call dista2(j(2),j(3),dis,b(1),b(2),b(3))
|
|
|
|
poem=rla*rlb
|
|
tel=a(1)*b(1)+a(2)*b(2)+a(3)*b(3)
|
|
arg=tel/poem
|
|
arg2=arg*arg
|
|
s1ma22=1.0-arg2
|
|
if (s1ma22.lt.1.0d-10) s1ma22=1.0d-10
|
|
s1ma2=sqrt(s1ma22)
|
|
if (arg.gt.1.0) arg=1.0
|
|
if (arg.lt.-1.0) arg=-1.0
|
|
hl=acos(arg)
|
|
h(i1)=hl
|
|
**********************************************************************
|
|
* *
|
|
* Calculate derivative valency angle to cartesian coordinates *
|
|
* *
|
|
**********************************************************************
|
|
if (j(1).eq.ib(la,2)) then
|
|
do k1=1,3
|
|
dradc(k1,1)=drdc(la,k1,1)
|
|
dradc(k1,2)=drdc(la,k1,2)
|
|
end do
|
|
else
|
|
do k1=1,3
|
|
dradc(k1,1)=drdc(la,k1,2)
|
|
dradc(k1,2)=drdc(la,k1,1)
|
|
end do
|
|
end if
|
|
if (j(2).eq.ib(lb,2)) then
|
|
do k1=1,3
|
|
drbdc(k1,2)=drdc(lb,k1,1)
|
|
drbdc(k1,3)=drdc(lb,k1,2)
|
|
end do
|
|
else
|
|
do k1=1,3
|
|
drbdc(k1,2)=drdc(lb,k1,2)
|
|
drbdc(k1,3)=drdc(lb,k1,1)
|
|
end do
|
|
end if
|
|
do k1=1,3
|
|
do k2=1,3
|
|
dndc(k1,k2)=rla*drbdc(k1,k2)+rlb*dradc(k1,k2)
|
|
dtdc(k1,k2)=a(k1)*dbdc(k2)+b(k1)*dadc(k2)
|
|
dargdc(k1,k2)=(dtdc(k1,k2)-arg*dndc(k1,k2))/poem
|
|
dhdc(i1,k1,k2)=-dargdc(k1,k2)/s1ma2
|
|
end do
|
|
end do
|
|
|
|
10 continue
|
|
|
|
return
|
|
end
|
|
**********************************************************************
|
|
**********************************************************************
|
|
|
|
subroutine boncor
|
|
|
|
**********************************************************************
|
|
#include "cbka.blk"
|
|
#include "cbkabo.blk"
|
|
#include "cbkc.blk"
|
|
#include "cbkbo.blk"
|
|
#include "cbkboncor.blk"
|
|
#include "cbkbosi.blk"
|
|
#include "cbkbopi.blk"
|
|
#include "cbkbopi2.blk"
|
|
#include "cbkconst.blk"
|
|
#include "cbkdbopi2ndc.blk"
|
|
#include "cbkdbopidc.blk"
|
|
#include "cbkdbopindc.blk"
|
|
#include "cbkff.blk"
|
|
#include "cbkia.blk"
|
|
#include "cbkidbo.blk"
|
|
#include "cbknubon2.blk"
|
|
#include "cbkrbo.blk"
|
|
#include "control.blk"
|
|
#include "small.blk"
|
|
#include "cbkdbodc.blk"
|
|
c$$$ if (ndebug.eq.1) then
|
|
c$$$C open (65,file='fort.65',status='unknown',access='append')
|
|
c$$$ write (65,*) 'In boncor'
|
|
c$$$ call timer(65)
|
|
c$$$ close (65)
|
|
c$$$ end if
|
|
**********************************************************************
|
|
* *
|
|
* Correction for overcoordination and 1-3 bond orders *
|
|
* *
|
|
**********************************************************************
|
|
**********************************************************************
|
|
* Description of variables used in this routine.
|
|
*
|
|
* ndebug: stored in cbka.blk; control-parameter
|
|
* i1: local do-loop counter
|
|
* i2: local do-loop counter
|
|
* k1: local do-loop counter
|
|
* k2: local do-loop counter
|
|
* nbon: stored in cbka.blk; number of bonds in system
|
|
* ibt: local integer; stores bond type
|
|
* ib(nbomax,3): stored in cbka.blk: bond distance identifiers
|
|
* j1: local integer; stores atom number 1st atom in bond
|
|
* j2: local integer; stores atom number 2nd atom in bond
|
|
* ovc(nbotym): stored in cbka.blk: force field parameter for
|
|
* overcoordination correction
|
|
* v13cor(nbotym): stored in cbka.blk: force field parameter for
|
|
* 1-3 bond order correction
|
|
* idbo1(nbomax): stored in cbka.blk; number of atoms in the
|
|
* derivative of the bond order
|
|
* idbo(nbomax,2*mbond): stored in cbka.blk; atom numbers of the
|
|
* atoms in the derivative of the bond order
|
|
* dbondc(nbomax,3,2*mbond): stored in cbka.blk; derivative of
|
|
* corrected total bond orders to cartesians
|
|
* dbosindc(nbomax,3,2*mbond): stored in cbka.blk; derivative of
|
|
* corrected sigma bond orders to cartesians
|
|
* dbopindc(nbomax,3,2*mbond): stored in cbka.blk; derivative of
|
|
* corrected pi bond orders to cartesians
|
|
* dbopi2ndc(nbomax,3,2*mbond): stored in cbka.blk; derivative of
|
|
* corrected double pi bond orders to cartesians
|
|
* dbodc(nbomax,3,2): stored in cbka.blk; derivative of
|
|
* uncorrected total bond orders to cartesians
|
|
* dbosidc(nbomax,3,2): stored in cbka.blk; derivative of
|
|
* uncorrected sigma bond orders to cartesians
|
|
* dbopidc(nbomax,3,2): stored in cbka.blk; derivative of
|
|
* uncorrected pi bond orders to cartesians
|
|
* dbopi2dc(nbomax,3,2): stored in cbka.blk; derivative of
|
|
* uncorrected double pi bond orders to cartesians
|
|
* boo: local variable; storage of uncorrected total bond order
|
|
* bo(nbomax): stored in cbka.blk; total bond order
|
|
* bopi(nbomax): stored in cbka.blk; pi bond order
|
|
* bopi2(nbomax): stored in cbka.blk; double pi bond order
|
|
* bopio: local variable; storage of uncorrected pi bond order
|
|
* bopi2o: local variable; storage of uncorrected double pi bond order
|
|
* iti: local integer; atom type first atom in bond
|
|
* itj: local integer; atom type second atom in bond
|
|
* ia(nat,mbond+3): stored in cbka.blk; connection table without bond
|
|
* order cutoff
|
|
* aboi: local variable: total bond order around atom i
|
|
* aboj: local variable: total bond order around atom j
|
|
* abo(nat): stored in cbka.blk; total bond order around atoms
|
|
* vp131: local variable; force field cross-term
|
|
* vp132: local variable; force field cross-term
|
|
* vp133: local variable; force field cross-term
|
|
* bo131(nsort): stored in cbka.blk; force field parameter for 1-3
|
|
* bond order correction
|
|
* bo132(nsort): stored in cbka.blk; force field parameter for 1-3
|
|
* bond order correction
|
|
* bo133(nsort): stored in cbka.blk; force field parameter for 1-3
|
|
* bond order correction
|
|
* corrtot:local variable; total correction on bond order
|
|
* dbodsboi1: local variable; derivative of bond order to sum of bond
|
|
* orders around atom i
|
|
* dbodsboj1: local variable; derivative of bond order to sum of bond
|
|
* orders around atom j
|
|
* ovi: local variable; overcoordination on atom i
|
|
* ovj: local variable; overcoordination on atom j
|
|
* aval(nat): stored in cbka.blk; nr. of valence electrons on atom
|
|
* exphu1: local variable; stores exponential
|
|
* exphu2: local variable; stores exponential
|
|
* exp11: local variable; stores exponential
|
|
* exp21: local variable; stores exponential
|
|
* vpar(npamax): stored in cbka.blk: general parameters
|
|
* exphu12: local variable; stores sum of exponential
|
|
* ovcor: local variable; temporary storage for BO/ovcor corr.
|
|
* huli: local variable; temporary storage for BO/ovcor corr.
|
|
* hulj: local variable; temporary storage for BO/ovcor corr.
|
|
* corr1: local variable; temporary storage for BO/ovcor corr.
|
|
* corr2: local variable; temporary storage for BO/ovcor corr.
|
|
* dbodsboi2: local variable; derivative of 1-3 BO correction to sum
|
|
* of bond orders around atom i
|
|
* dbodsboj2: local variable; derivative of 1-3 BO correction to sum
|
|
* of bond orders around atom i
|
|
* bocor1: local variable; 1-3 bond order correction
|
|
* bocor2: local variable; 1-3 bond order correction
|
|
* ovi2: local variable; overcoordination on atom i with reference to
|
|
* total number of electrons on atom i, including lone
|
|
* pairs
|
|
* ovj2: local variable; overcoordination on atom j with reference to
|
|
* total number of electrons on atom j, including lone
|
|
* pairs
|
|
* valf(nsort): stored in cbka.blk; total number of electrons on
|
|
* atom, including lone pairs
|
|
* cor1: local variable; temporary storage for BO/1-3 bond corr.
|
|
* cor2: local variable; temporary storage for BO/1-3 bond corr.
|
|
* exphu3: local variable; storage exponential
|
|
* exphu4: local variable; storage exponential
|
|
* corrtot2: local variable; square of corrtot
|
|
* dbodboo: local variable; derivative of corrected total bond order to
|
|
* uncorrected bond order
|
|
* dbopidbopio: local variable; derivative of corrected pi bond order
|
|
* to uncorrected pi bond order
|
|
* dbopidboo: local variable; derivative of corrected pi bond order
|
|
* to uncorrected total bond order
|
|
* dbopi2dbopi2o: local variable; derivative of corrected double pi bond order
|
|
* to uncorrected double pi bond order
|
|
* dbopi2dboo: local variable; derivative of corrected double pi bond order
|
|
* to uncorrected total bond order
|
|
* dbodsboit: local variable; derivative of total bond order to sum
|
|
* of bond orders around atom i
|
|
* dbodsbojt: local variable; derivative of total bond order to sum
|
|
* of bond orders around atom j
|
|
* vhui: local variable; temporary storage
|
|
* vhuj: local variable; temporary storage
|
|
* dbopidsboit: local variable; derivative of pi bond order to sum
|
|
* of bond orders around atom i
|
|
* dbopidsbojt: local variable; derivative of pi bond order to sum
|
|
* of bond orders around atom j
|
|
* dbopi2dsboit: local variable; derivative of pi bond order to sum
|
|
* of bond orders around atom i
|
|
* dbopi2dsbojt: local variable; derivative of pi bond order to sum
|
|
* of bond orders around atom j
|
|
* nco: local integer; counter for number of atoms in derivative
|
|
* ihl: local integer; helps to access right dbodc-term
|
|
* nubon2(nat,mbond): stored in cbka.blk; stored bond number as a
|
|
* function of atom number and connection number
|
|
* iob: local integer; atom number of second atom in bond
|
|
* ncubo: local integer; stores number of current bond
|
|
* na: stored in cbka.blk: number of atoms in system
|
|
* zero: stored in cbka.blk: value 0.00
|
|
*
|
|
**********************************************************************
|
|
do 10 i1=1,nbon
|
|
ibt=ib(i1,1)
|
|
j1=ib(i1,2)
|
|
j2=ib(i1,3)
|
|
if (ovc(ibt).lt.0.001.and.v13cor(ibt).lt.0.001) then
|
|
idbo1(i1)=2
|
|
idbo(i1,1)=j1
|
|
idbo(i1,2)=j2
|
|
do k1=1,3
|
|
dbondc(i1,k1,1)=dbodc(i1,k1,1)
|
|
dbondc(i1,k1,2)=dbodc(i1,k1,2)
|
|
dbosindc(i1,k1,1)=dbosidc(i1,k1,1)
|
|
dbosindc(i1,k1,2)=dbosidc(i1,k1,2)
|
|
dbopindc(i1,k1,1)=dbopidc(i1,k1,1)
|
|
dbopindc(i1,k1,2)=dbopidc(i1,k1,2)
|
|
dbopi2ndc(i1,k1,1)=dbopi2dc(i1,k1,1)
|
|
dbopi2ndc(i1,k1,2)=dbopi2dc(i1,k1,2)
|
|
end do
|
|
goto 10
|
|
end if
|
|
boo=bo(i1)
|
|
bopio=bopi(i1)
|
|
bopi2o=bopi2(i1)
|
|
iti=ia(j1,1)
|
|
itj=ia(j2,1)
|
|
aboi=abo(j1)
|
|
aboj=abo(j2)
|
|
vp131=sqrt(bo131(iti)*bo131(itj))
|
|
vp132=sqrt(bo132(iti)*bo132(itj))
|
|
vp133=sqrt(bo133(iti)*bo133(itj))
|
|
corrtot=1.0
|
|
dbodsboi1=zero
|
|
dbodsboj1=zero
|
|
if (ovc(ibt).gt.0.001) then
|
|
ovi=aboi-aval(iti)
|
|
ovj=aboj-aval(itj)
|
|
|
|
**********************************************************************
|
|
* *
|
|
* Correction for overcoordination *
|
|
* *
|
|
**********************************************************************
|
|
exphu1=exp(-vpar(2)*ovi)
|
|
exphu2=exp(-vpar(2)*ovj)
|
|
exp11=exp(-vpar(1)*ovi)
|
|
exp21=exp(-vpar(1)*ovj)
|
|
exphu12=(exphu1+exphu2)
|
|
ovcor=-(1.0/vpar(2))*log(0.50*exphu12)
|
|
* huli=((1.0/ovc(ibt))*aval(iti)+exp11+exp21)
|
|
* hulj=((1.0/ovc(ibt))*aval(itj)+exp11+exp21)
|
|
huli=aval(iti)+exp11+exp21
|
|
hulj=aval(itj)+exp11+exp21
|
|
corr1=huli/(huli+ovcor)
|
|
corr2=hulj/(hulj+ovcor)
|
|
corrtot=0.50*(corr1+corr2)
|
|
|
|
dbodsboi1=0.50*(-vpar(1)*exp11/(huli+ovcor)-
|
|
$(corr1/(huli+ovcor))*
|
|
$(-vpar(1)*exp11+exphu1/exphu12)-vpar(1)*exp11/(hulj+ovcor)-
|
|
$(corr2/(hulj+ovcor))*(-vpar(1)*exp11+exphu1/exphu12))
|
|
dbodsboj1=0.50*(-vpar(1)*exp21/(huli+ovcor)-
|
|
$(corr1/(huli+ovcor))*
|
|
$(-vpar(1)*exp21+exphu2/exphu12)-vpar(1)*exp21/(hulj+ovcor)-
|
|
$(corr2/(hulj+ovcor))*(-vpar(1)*exp21+exphu2/exphu12))
|
|
end if
|
|
**********************************************************************
|
|
* *
|
|
* Correction for 1-3 bond orders *
|
|
* *
|
|
**********************************************************************
|
|
dbodsboi2=zero
|
|
dbodsboj2=zero
|
|
bocor1=1.0
|
|
bocor2=1.0
|
|
if (v13cor(ibt).gt.0.001) then
|
|
ovi2=aboi-vval3(iti) !Modification for metal surfaces
|
|
ovj2=aboj-vval3(itj)
|
|
* ovi2=aboi-valf(iti)
|
|
* ovj2=aboj-valf(itj)
|
|
* ovi2=aboi-aval(iti)
|
|
* ovj2=aboj-aval(itj)
|
|
cor1=vp131*boo*boo-ovi2
|
|
cor2=vp131*boo*boo-ovj2
|
|
* exphu3=v13cor(ibt)*exp(-vp132*cor1+vp133)
|
|
* exphu4=v13cor(ibt)*exp(-vp132*cor2+vp133)
|
|
exphu3=exp(-vp132*cor1+vp133)
|
|
exphu4=exp(-vp132*cor2+vp133)
|
|
bocor1=1.0/(1.0+exphu3)
|
|
bocor2=1.0/(1.0+exphu4)
|
|
dbodsboi2=-bocor1*bocor1*bocor2*vp132*exphu3
|
|
dbodsboj2=-bocor1*bocor2*bocor2*vp132*exphu4
|
|
end if
|
|
|
|
bo(i1)=boo*corrtot*bocor1*bocor2
|
|
if (bo(i1).lt.1e-10) bo(i1)=zero
|
|
corrtot2=corrtot*corrtot
|
|
bopi(i1)=bopio*corrtot2*bocor1*bocor2
|
|
bopi2(i1)=bopi2o*corrtot2*bocor1*bocor2
|
|
if (bopi(i1).lt.1e-10) bopi(i1)=zero
|
|
if (bopi2(i1).lt.1e-10) bopi2(i1)=zero
|
|
|
|
dbodboo=corrtot*bocor1*bocor2+corrtot*
|
|
$bocor1*bocor1*bocor2*boo*vp132*vp131*2.0*boo*exphu3+
|
|
$corrtot*bocor1*bocor2*bocor2*boo*
|
|
$vp132*vp131*exphu4*2.0*boo
|
|
|
|
dbopidbopio=corrtot2*bocor1*bocor2
|
|
|
|
dbopidboo=corrtot2*
|
|
$bocor1*bocor1*bocor2*boo*vp132*vp131*2.0*bopio*exphu3+
|
|
$corrtot2*bocor1*bocor2*bocor2*boo*
|
|
$vp132*vp131*exphu4*2.0*bopio
|
|
|
|
dbopi2dbopi2o=corrtot2*bocor1*bocor2
|
|
|
|
dbopi2dboo=corrtot2*
|
|
$bocor1*bocor1*bocor2*boo*vp132*vp131*2.0*bopi2o*exphu3+
|
|
$corrtot2*bocor1*bocor2*bocor2*boo*
|
|
$vp132*vp131*exphu4*2.0*bopi2o
|
|
|
|
dbodsboit=boo*dbodsboi1*bocor1*bocor2+boo*corrtot*dbodsboi2
|
|
dbodsbojt=boo*dbodsboj1*bocor1*bocor2+boo*corrtot*dbodsboj2
|
|
|
|
vhui=2.0*corrtot*dbodsboi1*bocor1*bocor2+corrtot2*dbodsboi2
|
|
vhuj=2.0*corrtot*dbodsboj1*bocor1*bocor2+corrtot2*dbodsboj2
|
|
dbopidsboit=bopio*vhui
|
|
dbopidsbojt=bopio*vhuj
|
|
|
|
dbopi2dsboit=bopi2o*vhui
|
|
dbopi2dsbojt=bopi2o*vhuj
|
|
|
|
**********************************************************************
|
|
* *
|
|
* Calculate bond order derivatives *
|
|
* *
|
|
**********************************************************************
|
|
idbo1(i1)=2+ia(j1,2)+ia(j2,2)
|
|
idbo(i1,1)=j1
|
|
idbo(i1,2)=j2
|
|
nco=0
|
|
do k1=1,3
|
|
dbondc(i1,k1,1)=dbodc(i1,k1,1)*dbodboo
|
|
dbondc(i1,k1,2)=dbodc(i1,k1,2)*dbodboo
|
|
* dbosindc(i1,k1,1)=dbosidc(i1,k1,1)*dbosidboo
|
|
* dbosindc(i1,k1,2)=dbosidc(i1,k1,2)*dbosidboo
|
|
dbopindc(i1,k1,1)=dbopidc(i1,k1,1)*dbopidbopio+
|
|
$dbodc(i1,k1,1)*dbopidboo
|
|
dbopindc(i1,k1,2)=dbopidc(i1,k1,2)*dbopidbopio+
|
|
$dbodc(i1,k1,2)*dbopidboo
|
|
dbopi2ndc(i1,k1,1)=dbopi2dc(i1,k1,1)*dbopi2dbopi2o+
|
|
$dbodc(i1,k1,1)*dbopi2dboo
|
|
dbopi2ndc(i1,k1,2)=dbopi2dc(i1,k1,2)*dbopi2dbopi2o+
|
|
$dbodc(i1,k1,2)*dbopi2dboo
|
|
end do
|
|
do i2=1,ia(j1,2)
|
|
ihl=0
|
|
iob=ia(j1,2+i2)
|
|
if (iob.lt.j1) ihl=1
|
|
ncubo=nubon2(j1,i2)
|
|
idbo(i1,2+nco+1)=iob
|
|
do k1=1,3
|
|
dbondc(i1,k1,1)=dbondc(i1,k1,1)+dbodc(ncubo,k1,1+ihl)*dbodsboit
|
|
dbondc(i1,k1,2+nco+1)=dbodc(ncubo,k1,2-ihl)*dbodsboit
|
|
|
|
* dbosindc(i1,k1,1)=dbosindc(i1,k1,1)+
|
|
* $dbodc(ncubo,k1,1+ihl)*dbosidsboit
|
|
* dbosindc(i1,k1,2+nco+1)=dbodc(ncubo,k1,2-ihl)*dbosidsboit
|
|
|
|
dbopindc(i1,k1,1)=dbopindc(i1,k1,1)+
|
|
$dbodc(ncubo,k1,1+ihl)*dbopidsboit
|
|
dbopindc(i1,k1,2+nco+1)=dbodc(ncubo,k1,2-ihl)*dbopidsboit
|
|
|
|
dbopi2ndc(i1,k1,1)=dbopi2ndc(i1,k1,1)+
|
|
$dbodc(ncubo,k1,1+ihl)*dbopi2dsboit
|
|
dbopi2ndc(i1,k1,2+nco+1)=dbodc(ncubo,k1,2-ihl)*dbopi2dsboit
|
|
|
|
end do
|
|
nco=nco+1
|
|
end do
|
|
do i2=1,ia(j2,2)
|
|
ihl=0
|
|
iob=ia(j2,2+i2)
|
|
if (iob.lt.j2) ihl=1
|
|
ncubo=nubon2(j2,i2)
|
|
idbo(i1,2+nco+1)=iob
|
|
do k1=1,3
|
|
|
|
dbondc(i1,k1,2)=dbondc(i1,k1,2)+dbodc(ncubo,k1,1+ihl)*dbodsbojt
|
|
dbondc(i1,k1,2+nco+1)=dbodc(ncubo,k1,2-ihl)*dbodsbojt
|
|
|
|
* dbosindc(i1,k1,2)=dbosindc(i1,k1,2)+
|
|
* $dbodc(ncubo,k1,1+ihl)*dbosidsbojt
|
|
* dbosindc(i1,k1,2+nco+1)=dbodc(ncubo,k1,2-ihl)*dbosidsbojt
|
|
|
|
dbopindc(i1,k1,2)=dbopindc(i1,k1,2)+
|
|
$dbodc(ncubo,k1,1+ihl)*dbopidsbojt
|
|
dbopindc(i1,k1,2+nco+1)=dbodc(ncubo,k1,2-ihl)*dbopidsbojt
|
|
|
|
dbopi2ndc(i1,k1,2)=dbopi2ndc(i1,k1,2)+
|
|
$dbodc(ncubo,k1,1+ihl)*dbopi2dsbojt
|
|
dbopi2ndc(i1,k1,2+nco+1)=dbodc(ncubo,k1,2-ihl)*dbopi2dsbojt
|
|
|
|
end do
|
|
nco=nco+1
|
|
end do
|
|
|
|
10 continue
|
|
|
|
do i1=1,na
|
|
abo(i1)=zero
|
|
end do
|
|
* do i1=1,na
|
|
* do i2=1,ia(i1,2)
|
|
* iob=ia(i1,2+i2)
|
|
* ncubo=nubon2(i1,i2)
|
|
* abo(i1)=abo(i1)+bo(ncubo)
|
|
* end do
|
|
* end do
|
|
do i1=1,nbon
|
|
j1=ib(i1,2)
|
|
j2=ib(i1,3)
|
|
abo(j1)=abo(j1)+bo(i1)
|
|
if (j1.ne.j2) abo(j2)=abo(j2)+bo(i1)
|
|
end do
|
|
|
|
15 continue
|
|
return
|
|
end
|
|
**********************************************************************
|
|
**********************************************************************
|
|
|
|
subroutine lonpar
|
|
|
|
**********************************************************************
|
|
#include "cbka.blk"
|
|
#include "cbkabo.blk"
|
|
#include "cbkconst.blk"
|
|
#include "cbkc.blk"
|
|
#include "cbkd.blk"
|
|
#include "cbkdcell.blk"
|
|
#include "cbkenergies.blk"
|
|
#include "cbkff.blk"
|
|
#include "cbkia.blk"
|
|
#include "cbkidbo.blk"
|
|
#include "cbklonpar.blk"
|
|
#include "cbknubon2.blk"
|
|
#include "control.blk"
|
|
#include "small.blk"
|
|
dimension virial_tmp(3,3),virialsym(6)
|
|
c$$$ if (ndebug.eq.1) then
|
|
c$$$C open (65,file='fort.65',status='unknown',access='append')
|
|
c$$$ write (65,*) 'In lonpar'
|
|
c$$$ call timer(65)
|
|
c$$$ close (65)
|
|
c$$$ end if
|
|
**********************************************************************
|
|
* *
|
|
* Calculate lone pair energy and first derivatives *
|
|
* *
|
|
**********************************************************************
|
|
elp=zero
|
|
do i1=1,na
|
|
**********************************************************************
|
|
* *
|
|
* Determine number of lone pairs on atoms
|
|
* *
|
|
**********************************************************************
|
|
ity=ia(i1,1)
|
|
voptlp=0.50*(stlp(ity)-aval(ity))
|
|
vlp(i1)=zero
|
|
vund=abo(i1)-stlp(ity)
|
|
vlph=2.0*int(vund/2.0)
|
|
vlpex=vund-vlph
|
|
vp16h=vpar(16)-1.0
|
|
|
|
expvlp=exp(-vpar(16)*(2.0+vlpex)*(2.0+vlpex))
|
|
dvlpdsbo(i1)=-vpar(16)*2.0*(2.0+vlpex)*expvlp
|
|
vlp(i1)=expvlp-int(vund/2.0)
|
|
* expvlp=exp(-vpar(16)*(2.0+vlpex))
|
|
* dvlpdsbo(i1)=-vpar(16)*expvlp
|
|
* expvlp=exp(-6.0*((-0.50*vlpex)**vpar(16)))
|
|
* vlp(i1)=(1.0-expvlp)-int(vund/2.0)
|
|
* dvlpdsbo(i1)=-0.5*6.0*vpar(16)*((-0.5*vlpex)**vp16h)*
|
|
* $expvlp
|
|
**********************************************************************
|
|
* *
|
|
* Calculate lone pair energy *
|
|
* *
|
|
**********************************************************************
|
|
if (i1 .le. na_local) then
|
|
|
|
diffvlp=voptlp-vlp(i1)
|
|
exphu1=exp(-75.0*diffvlp)
|
|
hulp1=1.0/(1.0+exphu1)
|
|
elph=vlp1(ity)*diffvlp*hulp1
|
|
* elph=vlp1(ity)*diffvlp
|
|
delpdvlp=-vlp1(ity)*hulp1-vlp1(ity)*diffvlp*hulp1*hulp1*
|
|
$75.0*exphu1
|
|
|
|
elp=elp+elph
|
|
estrain(i1)=estrain(i1)+elph !atom energy
|
|
|
|
delpdsbo=delpdvlp*dvlpdsbo(i1)
|
|
**********************************************************************
|
|
* *
|
|
* Calculate first derivative of lone pair energy to *
|
|
* cartesian coordinates *
|
|
* *
|
|
**********************************************************************
|
|
do i3=1,ia(i1,2)
|
|
iob=ia(i1,2+i3)
|
|
ncubo=nubon2(i1,i3)
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1=1,3
|
|
do k2=1,3
|
|
virial_tmp(k1,k2) = 0.0
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
do i4=1,idbo1(ncubo)
|
|
ihu=idbo(ncubo,i4)
|
|
do k1=1,3
|
|
ftmp = delpdsbo*dbondc(ncubo,k1,i4)
|
|
d(k1,ihu)=d(k1,ihu)+ftmp
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1p=1,3
|
|
virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p)
|
|
end do
|
|
endif
|
|
|
|
end do
|
|
end do
|
|
if (Lvirial.eq.1) then
|
|
virialsym(1) = virial_tmp(1,1)
|
|
virialsym(2) = virial_tmp(2,2)
|
|
virialsym(3) = virial_tmp(3,3)
|
|
virialsym(4) = virial_tmp(1,2)
|
|
virialsym(5) = virial_tmp(1,3)
|
|
virialsym(6) = virial_tmp(2,3)
|
|
do k1 = 1,6
|
|
virial(k1) = virial(k1) + virialsym(k1)
|
|
end do
|
|
|
|
if (Latomvirial.eq.1) then
|
|
frac = 1.0d0/idbo1(ncubo)
|
|
do k1 = 1,6
|
|
vtmp = virialsym(k1)*frac
|
|
do i4=1,idbo1(ncubo)
|
|
ihu=idbo(ncubo,i4)
|
|
atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp
|
|
end do
|
|
end do
|
|
endif
|
|
endif
|
|
|
|
end do
|
|
endif
|
|
|
|
end do
|
|
|
|
return
|
|
end
|
|
**********************************************************************
|
|
**********************************************************************
|
|
|
|
subroutine covbon
|
|
|
|
**********************************************************************
|
|
#include "cbka.blk"
|
|
#include "cbkc.blk"
|
|
#include "cbkabo.blk"
|
|
#include "cbkbo.blk"
|
|
#include "cbkbosi.blk"
|
|
#include "cbkbopi.blk"
|
|
#include "cbkbopi2.blk"
|
|
#include "cbkconst.blk"
|
|
#include "cbkcovbon.blk"
|
|
#include "cbkd.blk"
|
|
#include "cbkdbopi2ndc.blk"
|
|
#include "cbkdbopindc.blk"
|
|
#include "cbkdcell.blk"
|
|
#include "cbkenergies.blk"
|
|
#include "cbkff.blk"
|
|
#include "cbkia.blk"
|
|
#include "cbkidbo.blk"
|
|
#include "cbknubon2.blk"
|
|
#include "cbkqa.blk"
|
|
#include "cbkrbo.blk"
|
|
#include "control.blk"
|
|
#include "small.blk"
|
|
dimension virial_tmp(3,3),virialsym(6)
|
|
**********************************************************************
|
|
* *
|
|
* Calculate bond energy and first derivatives *
|
|
* *
|
|
**********************************************************************
|
|
c$$$ if (ndebug.eq.1) then
|
|
c$$$C open (65,file='fort.65',status='unknown',access='append')
|
|
c$$$ write (65,*) 'In covbon'
|
|
c$$$ call timer(65)
|
|
c$$$ close (65)
|
|
c$$$ end if
|
|
eb=0.0d0
|
|
if (nbon.eq.0) return
|
|
**********************************************************************
|
|
* *
|
|
* Calculate bond energies *
|
|
* *
|
|
**********************************************************************
|
|
c$$$ if (ndebug.eq.1) then
|
|
c$$$C open (65,file='fort.65',status='unknown',access='append')
|
|
c$$$ write(65,*) 'Bond forces'
|
|
c$$$ write(65,*) 'nbon = ',nbon
|
|
c$$$ endif
|
|
|
|
do 20 i1=1,nbon
|
|
|
|
boa=bo(i1)
|
|
* if (boa.lt.cutof2) goto 20
|
|
j1=ib(i1,2)
|
|
j2=ib(i1,3)
|
|
|
|
c Only compute interaction if both atoms
|
|
c are local or else flip a coin
|
|
if (j1 .gt. na_local) go to 20
|
|
if (j2 .gt. na_local) then
|
|
if (itag(j1) .lt. itag(j2)) go to 20
|
|
if (itag(j1) .eq. itag(j2)) then
|
|
if(c(j1,3) .gt. c(j2,3)) go to 20
|
|
if(c(j1,3) .eq. c(j2,3) .and.
|
|
$ c(j1,2) .gt. c(j2,2)) go to 20
|
|
if(c(j1,3) .eq. c(j2,3) .and.
|
|
$ c(j1,2) .eq. c(j2,2) .and.
|
|
$ c(j1,1) .gt. c(j2,1)) go to 20
|
|
endif
|
|
endif
|
|
vsymm=1.0
|
|
if (j1.eq.j2) vsymm=0.5
|
|
|
|
bopia=bopi(i1)
|
|
bopi2a=bopi2(i1)
|
|
bosia=boa-bopia-bopi2a
|
|
if (bosia.lt.zero) bosia=zero
|
|
it1=ia(j1,1)
|
|
it2=ia(j2,1)
|
|
ibt=ib(i1,1)
|
|
de1h=vsymm*de1(ibt)
|
|
de2h=vsymm*de2(ibt)
|
|
de3h=vsymm*de3(ibt)
|
|
|
|
bopo1=bosia**psp(ibt)
|
|
exphu1=exp(psi(ibt)*(1.0-bopo1))
|
|
ebh=-de1h*bosia*exphu1-de2h*bopia-de3h*bopi2a
|
|
|
|
debdbo=-de1h*exphu1+de1h*exphu1*psp(ibt)*psi(ibt)*bopo1
|
|
debdbopi=-de2h
|
|
debdbopi2=-de3h
|
|
|
|
eb=eb+ebh
|
|
estrain(j1)=estrain(j1)+0.50*ebh !1st atom energy
|
|
estrain(j2)=estrain(j2)+0.50*ebh !2nd atom energy
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1=1,3
|
|
do k2=1,3
|
|
virial_tmp(k1,k2) = 0.0
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
do i2=1,idbo1(i1)
|
|
ihu=idbo(i1,i2)
|
|
do k1=1,3
|
|
ftmp = debdbo*(dbondc(i1,k1,i2)-dbopindc(i1,k1,i2)-
|
|
$dbopi2ndc(i1,k1,i2))+
|
|
$debdbopi*dbopindc(i1,k1,i2)+
|
|
$debdbopi2*dbopi2ndc(i1,k1,i2)
|
|
d(k1,ihu)=d(k1,ihu)+ftmp
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1p=1,3
|
|
virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p)
|
|
end do
|
|
endif
|
|
|
|
end do
|
|
end do
|
|
if (Lvirial.eq.1) then
|
|
virialsym(1) = virial_tmp(1,1)
|
|
virialsym(2) = virial_tmp(2,2)
|
|
virialsym(3) = virial_tmp(3,3)
|
|
virialsym(4) = virial_tmp(1,2)
|
|
virialsym(5) = virial_tmp(1,3)
|
|
virialsym(6) = virial_tmp(2,3)
|
|
do k1 = 1,6
|
|
virial(k1) = virial(k1) + virialsym(k1)
|
|
end do
|
|
|
|
if (Latomvirial.eq.1) then
|
|
frac = 1.0d0/idbo1(i1)
|
|
do k1 = 1,6
|
|
vtmp = virialsym(k1)*frac
|
|
do i2=1,idbo1(i1)
|
|
ihu=idbo(i1,i2)
|
|
atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
endif
|
|
|
|
**********************************************************************
|
|
* *
|
|
* Stabilisation terminal triple bond in CO *
|
|
* *
|
|
**********************************************************************
|
|
if (boa.lt.1.00) goto 20
|
|
* Stabilization for all triple bonds (not just for CO) in ReaxFF combustion FF
|
|
if (ltripstaball.eq.1 .or.
|
|
$ (qa(j1).eq.'C '.and.qa(j2).eq.'O ').or.
|
|
$ (qa(j1).eq.'O '.and.qa(j2).eq.'C ')) then
|
|
|
|
ba=(boa-2.50)*(boa-2.50)
|
|
exphu=exp(-vpar(8)*ba)
|
|
oboa=abo(j1)-boa
|
|
obob=abo(j2)-boa
|
|
exphua1=exp(-vpar(4)*oboa)
|
|
exphub1=exp(-vpar(4)*obob)
|
|
ovoab=abo(j1)-aval(it1)+abo(j2)-aval(it2)
|
|
exphuov=exp(vpar(5)*ovoab)
|
|
hulpov=1.0/(1.0+25.0*exphuov)
|
|
|
|
estriph=vpar(11)*exphu*hulpov*(exphua1+exphub1)
|
|
|
|
eb=eb+estriph
|
|
estrain(j1)=estrain(j1)+0.50*estriph !1st atom energy
|
|
estrain(j2)=estrain(j2)+0.50*estriph !2nd atom energy
|
|
|
|
decobdbo=vpar(4)*vpar(11)*exphu*hulpov*(exphua1+exphub1)
|
|
$-2.0*vpar(11)*vpar(8)*(boa-2.50)*hulpov*exphu*
|
|
$(exphua1+exphub1)
|
|
decobdboua=-25.0*vpar(5)*vpar(11)*exphu*exphuov*hulpov*hulpov*
|
|
$(exphua1+exphub1)-vpar(11)*exphu*vpar(4)*hulpov*exphua1
|
|
decobdboub=-25.0*vpar(5)*vpar(11)*exphu*exphuov*hulpov*hulpov*
|
|
$(exphua1+exphub1)-vpar(11)*exphu*vpar(4)*hulpov*exphub1
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1=1,3
|
|
do k2=1,3
|
|
virial_tmp(k1,k2) = 0.0
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
do i2=1,idbo1(i1)
|
|
ihu=idbo(i1,i2)
|
|
do k1=1,3
|
|
ftmp = decobdbo*dbondc(i1,k1,i2)
|
|
d(k1,ihu)=d(k1,ihu)+ftmp
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1p=1,3
|
|
virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p)
|
|
end do
|
|
endif
|
|
|
|
end do
|
|
end do
|
|
if (Lvirial.eq.1) then
|
|
virialsym(1) = virial_tmp(1,1)
|
|
virialsym(2) = virial_tmp(2,2)
|
|
virialsym(3) = virial_tmp(3,3)
|
|
virialsym(4) = virial_tmp(1,2)
|
|
virialsym(5) = virial_tmp(1,3)
|
|
virialsym(6) = virial_tmp(2,3)
|
|
do k1 = 1,6
|
|
virial(k1) = virial(k1) + virialsym(k1)
|
|
end do
|
|
|
|
if (Latomvirial.eq.1) then
|
|
frac = 1.0d0/idbo1(i1)
|
|
do k1 = 1,6
|
|
vtmp = virialsym(k1)*frac
|
|
do i2=1,idbo1(i1)
|
|
ihu=idbo(i1,i2)
|
|
atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
endif
|
|
|
|
do i3=1,ia(j1,2)
|
|
iob=ia(j1,2+i3)
|
|
ncubo=nubon2(j1,i3)
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1=1,3
|
|
do k2=1,3
|
|
virial_tmp(k1,k2) = 0.0
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
do i4=1,idbo1(ncubo)
|
|
ihu=idbo(ncubo,i4)
|
|
do k1=1,3
|
|
ftmp = decobdboua*dbondc(ncubo,k1,i4)
|
|
d(k1,ihu)=d(k1,ihu)+ftmp
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1p=1,3
|
|
virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p)
|
|
end do
|
|
endif
|
|
|
|
end do
|
|
end do
|
|
if (Lvirial.eq.1) then
|
|
virialsym(1) = virial_tmp(1,1)
|
|
virialsym(2) = virial_tmp(2,2)
|
|
virialsym(3) = virial_tmp(3,3)
|
|
virialsym(4) = virial_tmp(1,2)
|
|
virialsym(5) = virial_tmp(1,3)
|
|
virialsym(6) = virial_tmp(2,3)
|
|
do k1 = 1,6
|
|
virial(k1) = virial(k1) + virialsym(k1)
|
|
end do
|
|
|
|
if (Latomvirial.eq.1) then
|
|
frac = 1.0d0/idbo1(ncubo)
|
|
do k1 = 1,6
|
|
vtmp = virialsym(k1)*frac
|
|
do i4=1,idbo1(ncubo)
|
|
ihu=idbo(ncubo,i4)
|
|
atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
endif
|
|
|
|
end do
|
|
|
|
do i3=1,ia(j2,2)
|
|
iob=ia(j2,2+i3)
|
|
ncubo=nubon2(j2,i3)
|
|
if (Lvirial.eq.1) then
|
|
do k1=1,3
|
|
do k2=1,3
|
|
virial_tmp(k1,k2) = 0.0
|
|
end do
|
|
end do
|
|
endif
|
|
do i4=1,idbo1(ncubo)
|
|
ihu=idbo(ncubo,i4)
|
|
do k1=1,3
|
|
ftmp = decobdboub*dbondc(ncubo,k1,i4)
|
|
d(k1,ihu)=d(k1,ihu)+ftmp
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1p=1,3
|
|
virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p)
|
|
end do
|
|
endif
|
|
|
|
end do
|
|
end do
|
|
if (Lvirial.eq.1) then
|
|
virialsym(1) = virial_tmp(1,1)
|
|
virialsym(2) = virial_tmp(2,2)
|
|
virialsym(3) = virial_tmp(3,3)
|
|
virialsym(4) = virial_tmp(1,2)
|
|
virialsym(5) = virial_tmp(1,3)
|
|
virialsym(6) = virial_tmp(2,3)
|
|
do k1 = 1,6
|
|
virial(k1) = virial(k1) + virialsym(k1)
|
|
end do
|
|
|
|
if (Latomvirial.eq.1) then
|
|
frac = 1.0d0/idbo1(ncubo)
|
|
do k1 = 1,6
|
|
vtmp = virialsym(k1)*frac
|
|
do i4=1,idbo1(ncubo)
|
|
ihu=idbo(ncubo,i4)
|
|
atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
endif
|
|
|
|
end do
|
|
|
|
endif
|
|
|
|
20 continue
|
|
|
|
return
|
|
end
|
|
**********************************************************************
|
|
**********************************************************************
|
|
|
|
subroutine ovcor
|
|
|
|
**********************************************************************
|
|
#include "cbka.blk"
|
|
#include "cbkc.blk"
|
|
#include "cbkabo.blk"
|
|
#include "cbkbo.blk"
|
|
#include "cbkbopi.blk"
|
|
#include "cbkbopi2.blk"
|
|
#include "cbkconst.blk"
|
|
#include "cbkd.blk"
|
|
#include "cbkdbopi2ndc.blk"
|
|
#include "cbkdbopindc.blk"
|
|
#include "cbkdcell.blk"
|
|
#include "cbkenergies.blk"
|
|
#include "cbkff.blk"
|
|
#include "cbkia.blk"
|
|
#include "cbkidbo.blk"
|
|
#include "cbklonpar.blk"
|
|
#include "cbknubon2.blk"
|
|
#include "cbkrbo.blk"
|
|
#include "control.blk"
|
|
#include "small.blk"
|
|
**********************************************************************
|
|
* *
|
|
* Calculate atom energy *
|
|
* Correction for over- and undercoordinated atoms *
|
|
* *
|
|
**********************************************************************
|
|
dimension vlptemp(nat)
|
|
dimension virial_tmp(3,3),virialsym(6)
|
|
c$$$ if (ndebug.eq.1) then
|
|
c$$$C open (65,file='fort.65',status='unknown',access='append')
|
|
c$$$ write (65,*) 'In ovcor'
|
|
c$$$ call timer(65)
|
|
c$$$ close (65)
|
|
c$$$ end if
|
|
do i1=1,na
|
|
ity1=ia(i1,1)
|
|
vlptemp(i1)=vlp(i1)
|
|
if (amas(ity1).gt.21.0) vlptemp(i1)=0.50*(stlp(ity1)-aval(ity1)) !Only for 1st-row elements
|
|
end do
|
|
25 ea=zero
|
|
eaot=zero
|
|
eaut=zero
|
|
epen=0.0
|
|
|
|
do 30 i1=1,na_local
|
|
ity1=ia(i1,1)
|
|
dfvl=1.0
|
|
if (amas(ity1).gt.21.0) dfvl=0.0 !Only for 1st-row elements
|
|
**********************************************************************
|
|
* *
|
|
* Calculate overcoordination energy *
|
|
* Valency is corrected for lone pairs *
|
|
* *
|
|
**********************************************************************
|
|
|
|
voptlp=0.50*(stlp(ity1)-aval(ity1))
|
|
diffvlph=dfvl*(voptlp-vlptemp(i1))
|
|
**********************************************************************
|
|
* *
|
|
* Determine coordination neighboring atoms *
|
|
* *
|
|
**********************************************************************
|
|
sumov=0.0
|
|
sumov2=0.0
|
|
do i3=1,ia(i1,2)
|
|
iat2=ia(i1,2+i3)
|
|
ity2=ia(iat2,1)
|
|
ncubo=nubon2(i1,i3)
|
|
if (bo(ncubo).gt.0.0) then
|
|
ibt=ib(ncubo,1)
|
|
voptlp2=0.50*(stlp(ity2)-aval(ity2))
|
|
diffvlp2=dfvl*(voptlp2-vlptemp(iat2))
|
|
sumov=sumov+(bopi(ncubo)+bopi2(ncubo))*
|
|
$(abo(iat2)-aval(ity2)-diffvlp2)
|
|
sumov2=sumov2+vover(ibt)*de1(ibt)*bo(ncubo)
|
|
endif
|
|
end do
|
|
|
|
exphu1=exp(vpar(32)*sumov)
|
|
vho=1.0/(1.0+vpar(33)*exphu1)
|
|
diffvlp=diffvlph*vho
|
|
|
|
vov1=abo(i1)-aval(ity1)-diffvlp
|
|
dvov1dsumov=diffvlph*vpar(32)*vpar(33)*vho*vho*exphu1
|
|
exphuo=exp(vovun(ity1)*vov1)
|
|
hulpo=1.0/(1.0+exphuo)
|
|
|
|
hulpp=(1.0/(vov1+aval(ity1)+1e-8))
|
|
|
|
eah=sumov2*hulpp*hulpo*vov1
|
|
deadvov1=-sumov2*hulpp*hulpp*vov1*hulpo+
|
|
$sumov2*hulpp*hulpo-sumov2*hulpp*vov1*vovun(ity1)*
|
|
$hulpo*hulpo*exphuo
|
|
|
|
ea=ea+eah
|
|
estrain(i1)=estrain(i1)+eah !atom energy
|
|
**********************************************************************
|
|
* *
|
|
* Calculate first derivative of overcoordination energy to *
|
|
* cartesian coordinates *
|
|
* *
|
|
**********************************************************************
|
|
do i3=1,ia(i1,2)
|
|
iob=ia(i1,2+i3)
|
|
ncubo=nubon2(i1,i3)
|
|
if (bo(ncubo).gt.0.0) then
|
|
ibt=ib(ncubo,1)
|
|
deadbo=vover(ibt)*de1(ibt)*hulpp*hulpo*vov1
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1=1,3
|
|
do k2=1,3
|
|
virial_tmp(k1,k2) = 0.0
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
do i4=1,idbo1(ncubo)
|
|
ihu=idbo(ncubo,i4)
|
|
do k1=1,3
|
|
ftmp = deadvov1*(1.0+dfvl*vho*dvlpdsbo(i1))*
|
|
$dbondc(ncubo,k1,i4)+deadbo*dbondc(ncubo,k1,i4)
|
|
d(k1,ihu)=d(k1,ihu)+ftmp
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1p=1,3
|
|
virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p)
|
|
end do
|
|
endif
|
|
|
|
end do
|
|
end do
|
|
if (Lvirial.eq.1) then
|
|
virialsym(1) = virial_tmp(1,1)
|
|
virialsym(2) = virial_tmp(2,2)
|
|
virialsym(3) = virial_tmp(3,3)
|
|
virialsym(4) = virial_tmp(1,2)
|
|
virialsym(5) = virial_tmp(1,3)
|
|
virialsym(6) = virial_tmp(2,3)
|
|
do k1 = 1,6
|
|
virial(k1) = virial(k1) + virialsym(k1)
|
|
end do
|
|
|
|
if (Latomvirial.eq.1) then
|
|
frac = 1.0d0/idbo1(ncubo)
|
|
do k1 = 1,6
|
|
vtmp = virialsym(k1)*frac
|
|
do i4=1,idbo1(ncubo)
|
|
ihu=idbo(ncubo,i4)
|
|
atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
endif
|
|
|
|
endif
|
|
end do
|
|
|
|
do i2=1,ia(i1,2)
|
|
|
|
iat2=ia(i1,2+i2)
|
|
ity2=ia(iat2,1)
|
|
nbosa=nubon2(i1,i2)
|
|
if (bo(nbosa).gt.0.0) then
|
|
deadvov2=deadvov1*dvov1dsumov*(bopi(nbosa)+bopi2(nbosa))
|
|
|
|
voptlp2=0.50*(stlp(ity2)-aval(ity2))
|
|
diffvlp2=dfvl*(voptlp2-vlptemp(iat2))
|
|
deadpibo=deadvov1*dvov1dsumov*(abo(iat2)-aval(ity2)-diffvlp2)
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1=1,3
|
|
do k2=1,3
|
|
virial_tmp(k1,k2) = 0.0
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
do i4=1,idbo1(nbosa)
|
|
ihu=idbo(nbosa,i4)
|
|
do k1=1,3
|
|
ftmp = deadpibo*(dbopindc(nbosa,k1,i4)+
|
|
$dbopi2ndc(nbosa,k1,i4))
|
|
d(k1,ihu)=d(k1,ihu)+ftmp
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1p=1,3
|
|
virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p)
|
|
end do
|
|
endif
|
|
|
|
end do
|
|
end do
|
|
if (Lvirial.eq.1) then
|
|
virialsym(1) = virial_tmp(1,1)
|
|
virialsym(2) = virial_tmp(2,2)
|
|
virialsym(3) = virial_tmp(3,3)
|
|
virialsym(4) = virial_tmp(1,2)
|
|
virialsym(5) = virial_tmp(1,3)
|
|
virialsym(6) = virial_tmp(2,3)
|
|
do k1 = 1,6
|
|
virial(k1) = virial(k1) + virialsym(k1)
|
|
end do
|
|
|
|
if (Latomvirial.eq.1) then
|
|
frac = 1.0d0/idbo1(nbosa)
|
|
do k1 = 1,6
|
|
vtmp = virialsym(k1)*frac
|
|
do i4=1,idbo1(nbosa)
|
|
ihu=idbo(nbosa,i4)
|
|
atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp
|
|
end do
|
|
end do
|
|
endif
|
|
endif
|
|
|
|
do i3=1,ia(iat2,2)
|
|
iob=ia(iat2,2+i3)
|
|
ncubo=nubon2(iat2,i3)
|
|
if (bo(ncubo).gt.0.0) then
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1=1,3
|
|
do k2=1,3
|
|
virial_tmp(k1,k2) = 0.0
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
do i4=1,idbo1(ncubo)
|
|
ihu=idbo(ncubo,i4)
|
|
do k1=1,3
|
|
ftmp = deadvov2*(1.0+dfvl*dvlpdsbo(iat2))*
|
|
$dbondc(ncubo,k1,i4)
|
|
d(k1,ihu)=d(k1,ihu)+ftmp
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1p=1,3
|
|
virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p)
|
|
end do
|
|
endif
|
|
|
|
end do
|
|
end do
|
|
if (Lvirial.eq.1) then
|
|
virialsym(1) = virial_tmp(1,1)
|
|
virialsym(2) = virial_tmp(2,2)
|
|
virialsym(3) = virial_tmp(3,3)
|
|
virialsym(4) = virial_tmp(1,2)
|
|
virialsym(5) = virial_tmp(1,3)
|
|
virialsym(6) = virial_tmp(2,3)
|
|
do k1 = 1,6
|
|
virial(k1) = virial(k1) + virialsym(k1)
|
|
end do
|
|
|
|
if (Latomvirial.eq.1) then
|
|
frac = 1.0d0/idbo1(ncubo)
|
|
do k1 = 1,6
|
|
vtmp = virialsym(k1)*frac
|
|
do i4=1,idbo1(ncubo)
|
|
ihu=idbo(ncubo,i4)
|
|
atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
endif
|
|
|
|
endif
|
|
end do
|
|
|
|
endif
|
|
|
|
end do
|
|
|
|
**********************************************************************
|
|
* *
|
|
* Calculate undercoordination energy *
|
|
* *
|
|
**********************************************************************
|
|
if (valp1(ity1).lt.zero) goto 30 !skip undercoordination
|
|
exphu2=exp(vpar(10)*sumov)
|
|
vuhu1=1.0+vpar(9)*exphu2
|
|
hulpu2=1.0/vuhu1
|
|
|
|
exphu3=-exp(vpar(7)*vov1)
|
|
hulpu3=-(1.0+exphu3)
|
|
|
|
dise2=valp1(ity1)
|
|
exphuu=exp(-vovun(ity1)*vov1)
|
|
hulpu=1.0/(1.0+exphuu)
|
|
eahu=dise2*hulpu*hulpu2*hulpu3
|
|
deaudvov1=dise2*hulpu2*vovun(ity1)*hulpu*hulpu*exphuu*hulpu3-
|
|
$dise2*hulpu*hulpu2*vpar(7)*exphu3
|
|
|
|
ea=ea+eahu
|
|
estrain(i1)=estrain(i1)+eahu !atom energy
|
|
|
|
deaudsumov=-dise2*hulpu*vpar(9)*vpar(10)*hulpu3*exphu2*
|
|
$hulpu2*hulpu2
|
|
|
|
**********************************************************************
|
|
* *
|
|
* Calculate first derivative of atom energy to cartesian *
|
|
* coordinates *
|
|
* *
|
|
**********************************************************************
|
|
|
|
do i3=1,ia(i1,2)
|
|
iob=ia(i1,2+i3)
|
|
ncubo=nubon2(i1,i3)
|
|
if (bo(ncubo).gt.0.0) then
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1=1,3
|
|
do k2=1,3
|
|
virial_tmp(k1,k2) = 0.0
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
do i4=1,idbo1(ncubo)
|
|
ihu=idbo(ncubo,i4)
|
|
do k1=1,3
|
|
ftmp = deaudvov1*(1.0+dfvl*vho*dvlpdsbo(i1))*
|
|
$dbondc(ncubo,k1,i4)
|
|
d(k1,ihu)=d(k1,ihu)+ftmp
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1p=1,3
|
|
virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p)
|
|
end do
|
|
endif
|
|
|
|
end do
|
|
end do
|
|
if (Lvirial.eq.1) then
|
|
virialsym(1) = virial_tmp(1,1)
|
|
virialsym(2) = virial_tmp(2,2)
|
|
virialsym(3) = virial_tmp(3,3)
|
|
virialsym(4) = virial_tmp(1,2)
|
|
virialsym(5) = virial_tmp(1,3)
|
|
virialsym(6) = virial_tmp(2,3)
|
|
do k1 = 1,6
|
|
virial(k1) = virial(k1) + virialsym(k1)
|
|
end do
|
|
|
|
if (Latomvirial.eq.1) then
|
|
frac = 1.0d0/idbo1(ncubo)
|
|
do k1 = 1,6
|
|
vtmp = virialsym(k1)*frac
|
|
do i4=1,idbo1(ncubo)
|
|
ihu=idbo(ncubo,i4)
|
|
atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
endif
|
|
|
|
endif
|
|
end do
|
|
|
|
do i2=1,ia(i1,2)
|
|
|
|
iat2=ia(i1,2+i2)
|
|
ity2=ia(iat2,1)
|
|
nbosa=nubon2(i1,i2)
|
|
if (bo(nbosa).gt.0.0) then
|
|
deadvov2=(deaudsumov+dvov1dsumov*deaudvov1)*
|
|
$(bopi(nbosa)+bopi2(nbosa))
|
|
|
|
voptlp2=0.50*(stlp(ity2)-aval(ity2))
|
|
diffvlp2=dfvl*(voptlp2-vlptemp(iat2))
|
|
deadpibo1=(dvov1dsumov*deaudvov1+deaudsumov)*
|
|
$(abo(iat2)-aval(ity2)-diffvlp2)
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1=1,3
|
|
do k2=1,3
|
|
virial_tmp(k1,k2) = 0.0
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
do i4=1,idbo1(nbosa)
|
|
ihu=idbo(nbosa,i4)
|
|
do k1=1,3
|
|
ftmp = deadpibo1*
|
|
$(dbopindc(nbosa,k1,i4)+dbopi2ndc(nbosa,k1,i4))
|
|
d(k1,ihu)=d(k1,ihu)+ftmp
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1p=1,3
|
|
virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p)
|
|
end do
|
|
endif
|
|
|
|
end do
|
|
end do
|
|
if (Lvirial.eq.1) then
|
|
virialsym(1) = virial_tmp(1,1)
|
|
virialsym(2) = virial_tmp(2,2)
|
|
virialsym(3) = virial_tmp(3,3)
|
|
virialsym(4) = virial_tmp(1,2)
|
|
virialsym(5) = virial_tmp(1,3)
|
|
virialsym(6) = virial_tmp(2,3)
|
|
do k1 = 1,6
|
|
virial(k1) = virial(k1) + virialsym(k1)
|
|
end do
|
|
|
|
if (Latomvirial.eq.1) then
|
|
frac = 1.0d0/idbo1(nbosa)
|
|
do k1 = 1,6
|
|
vtmp = virialsym(k1)*frac
|
|
do i4=1,idbo1(nbosa)
|
|
ihu=idbo(nbosa,i4)
|
|
atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
endif
|
|
|
|
do i3=1,ia(iat2,2)
|
|
iob=ia(iat2,2+i3)
|
|
ncubo=nubon2(iat2,i3)
|
|
if (bo(ncubo).gt.0.0) then
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1=1,3
|
|
do k2=1,3
|
|
virial_tmp(k1,k2) = 0.0
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
do i4=1,idbo1(ncubo)
|
|
ihu=idbo(ncubo,i4)
|
|
do k1=1,3
|
|
ftmp = deadvov2*(1.0+dfvl*dvlpdsbo(iat2))*
|
|
$dbondc(ncubo,k1,i4)
|
|
d(k1,ihu)=d(k1,ihu)+ftmp
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1p=1,3
|
|
virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p)
|
|
end do
|
|
endif
|
|
|
|
end do
|
|
end do
|
|
if (Lvirial.eq.1) then
|
|
virialsym(1) = virial_tmp(1,1)
|
|
virialsym(2) = virial_tmp(2,2)
|
|
virialsym(3) = virial_tmp(3,3)
|
|
virialsym(4) = virial_tmp(1,2)
|
|
virialsym(5) = virial_tmp(1,3)
|
|
virialsym(6) = virial_tmp(2,3)
|
|
do k1 = 1,6
|
|
virial(k1) = virial(k1) + virialsym(k1)
|
|
end do
|
|
|
|
if (Latomvirial.eq.1) then
|
|
frac = 1.0d0/idbo1(ncubo)
|
|
do k1 = 1,6
|
|
vtmp = virialsym(k1)*frac
|
|
do i4=1,idbo1(ncubo)
|
|
ihu=idbo(ncubo,i4)
|
|
atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
endif
|
|
|
|
endif
|
|
end do
|
|
|
|
endif
|
|
|
|
end do
|
|
|
|
|
|
30 continue
|
|
|
|
**********************************************************************
|
|
* *
|
|
* Calculate correction for C2 *
|
|
* *
|
|
**********************************************************************
|
|
if (abs(vpar(6)).gt.0.001) then
|
|
do 40 i1=1,na_local
|
|
ity1=ia(i1,1)
|
|
vov4=abo(i1)-aval(ity1)
|
|
|
|
do i2=1,ia(i1,2)
|
|
iat2=ia(i1,2+i2)
|
|
nbohu=nubon2(i1,i2)
|
|
if (bo(nbohu).gt.0.0) then
|
|
|
|
ibt=ib(nbohu,1)
|
|
elph=zero
|
|
deahu2dbo=zero
|
|
deahu2dsbo=zero
|
|
vov3=bo(nbohu)-vov4-0.040*(vov4**4)
|
|
if (vov3.gt.3.0) then
|
|
elph=vpar(6)*(vov3-3.0)*(vov3-3.0)
|
|
deahu2dbo=2.0*vpar(6)*(vov3-3.0)
|
|
deahu2dsbo=2.0*vpar(6)*(vov3-3.0)*(-1.0-
|
|
$0.16*(vov4**3))
|
|
end if
|
|
|
|
elp=elp+elph
|
|
estrain(i1)=estrain(i1)+elph !atom energy
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1=1,3
|
|
do k2=1,3
|
|
virial_tmp(k1,k2) = 0.0
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
do i3=1,idbo1(nbohu)
|
|
ihu=idbo(nbohu,i3)
|
|
do k1=1,3
|
|
ftmp = deahu2dbo*dbondc(nbohu,k1,i3)
|
|
d(k1,ihu)=d(k1,ihu)+ftmp
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1p=1,3
|
|
virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p)
|
|
end do
|
|
endif
|
|
|
|
end do
|
|
end do
|
|
if (Lvirial.eq.1) then
|
|
virialsym(1) = virial_tmp(1,1)
|
|
virialsym(2) = virial_tmp(2,2)
|
|
virialsym(3) = virial_tmp(3,3)
|
|
virialsym(4) = virial_tmp(1,2)
|
|
virialsym(5) = virial_tmp(1,3)
|
|
virialsym(6) = virial_tmp(2,3)
|
|
do k1 = 1,6
|
|
virial(k1) = virial(k1) + virialsym(k1)
|
|
end do
|
|
|
|
if (Latomvirial.eq.1) then
|
|
frac = 1.0d0/idbo1(nbohu)
|
|
do k1 = 1,6
|
|
vtmp = virialsym(k1)*frac
|
|
do i3=1,idbo1(nbohu)
|
|
ihu=idbo(nbohu,i3)
|
|
atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
endif
|
|
|
|
do i3=1,ia(i1,2)
|
|
iob=ia(i1,2+i3)
|
|
ncubo=nubon2(i1,i3)
|
|
if (bo(ncubo).gt.0.0) then
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1=1,3
|
|
do k2=1,3
|
|
virial_tmp(k1,k2) = 0.0
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
do i4=1,idbo1(ncubo)
|
|
ihu=idbo(ncubo,i4)
|
|
do k1=1,3
|
|
ftmp = deahu2dsbo*dbondc(ncubo,k1,i4)
|
|
d(k1,ihu)=d(k1,ihu)+ftmp
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1p=1,3
|
|
virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p)
|
|
end do
|
|
endif
|
|
|
|
end do
|
|
end do
|
|
if (Lvirial.eq.1) then
|
|
virialsym(1) = virial_tmp(1,1)
|
|
virialsym(2) = virial_tmp(2,2)
|
|
virialsym(3) = virial_tmp(3,3)
|
|
virialsym(4) = virial_tmp(1,2)
|
|
virialsym(5) = virial_tmp(1,3)
|
|
virialsym(6) = virial_tmp(2,3)
|
|
do k1 = 1,6
|
|
virial(k1) = virial(k1) + virialsym(k1)
|
|
end do
|
|
|
|
if (Latomvirial.eq.1) then
|
|
frac = 1.0d0/idbo1(ncubo)
|
|
do k1 = 1,6
|
|
vtmp = virialsym(k1)*frac
|
|
do i4=1,idbo1(ncubo)
|
|
ihu=idbo(ncubo,i4)
|
|
atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
endif
|
|
|
|
end if
|
|
end do
|
|
|
|
end if
|
|
end do
|
|
|
|
40 continue
|
|
end if
|
|
|
|
return
|
|
end
|
|
**********************************************************************
|
|
**********************************************************************
|
|
|
|
subroutine molen
|
|
|
|
**********************************************************************
|
|
#include "cbka.blk"
|
|
#include "cbkbo.blk"
|
|
#include "cbkconst.blk"
|
|
#include "cbkc.blk"
|
|
#include "cbkd.blk"
|
|
#include "cbkenergies.blk"
|
|
#include "cbkff.blk"
|
|
#include "cbkia.blk"
|
|
#include "cbkidbo.blk"
|
|
#include "cbknmolat.blk"
|
|
#include "cbknubon2.blk"
|
|
#include "control.blk"
|
|
#include "small.blk"
|
|
dimension virial_tmp(3,3),virialsym(6)
|
|
**********************************************************************
|
|
* *
|
|
* Calculate molecular energy and first derivatives *
|
|
* Only used to prevent creating virtual electrons *
|
|
* *
|
|
**********************************************************************
|
|
c$$$ if (ndebug.eq.1) then
|
|
c$$$C open (65,file='fort.65',status='unknown',access='append')
|
|
c$$$ write (65,*) 'In molen'
|
|
c$$$ call timer(65)
|
|
c$$$ close (65)
|
|
c$$$ end if
|
|
emol=zero
|
|
return
|
|
do i1=1,nmolo
|
|
|
|
enelm=0.0
|
|
do i2=1,na
|
|
if (ia(i2,3+mbond).eq.i1) then
|
|
it1=ia(i2,1)
|
|
enelm=enelm+aval(it1)
|
|
end if
|
|
end do
|
|
|
|
na1m=nmolat(i1,1)
|
|
|
|
enelm=2*int(enelm*0.50)
|
|
* enelm=elmol(i1)
|
|
bomsum=zero
|
|
do i2=1,na1m
|
|
ihu=nmolat(i1,i2+1)
|
|
do i3=1,ia(ihu,2)
|
|
ihu2=nubon2(ihu,i3)
|
|
bomsum=bomsum+bo(ihu2)
|
|
end do
|
|
end do
|
|
diff=(bomsum-enelm)
|
|
exphu=exp(-vpar(37)*diff)
|
|
exphu2=1.0/(1.0+15.0*exphu)
|
|
emolh=zero
|
|
demoldsbo=zero
|
|
emolh=vpar(38)*exphu2
|
|
emol=emol+emolh
|
|
demoldsbo=vpar(38)*vpar(37)*15.0*exphu2*exphu2*exphu
|
|
|
|
do i2=1,na1m
|
|
ihu1=nmolat(i1,i2+1)
|
|
do i3=1,ia(ihu1,2)
|
|
iob=ia(ihu1,2+i3)
|
|
ncubo=nubon2(ihu1,i3)
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1=1,3
|
|
do k2=1,3
|
|
virial_tmp(k1,k2) = 0.0
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
do i4=1,idbo1(ncubo)
|
|
ihu=idbo(ncubo,i4)
|
|
do k1=1,3
|
|
ftmp = demoldsbo*dbondc(ncubo,k1,i4)
|
|
d(k1,ihu)=d(k1,ihu)+ftmp
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1p=1,3
|
|
virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p)
|
|
end do
|
|
endif
|
|
|
|
end do
|
|
end do
|
|
if (Lvirial.eq.1) then
|
|
virialsym(1) = virial_tmp(1,1)
|
|
virialsym(2) = virial_tmp(2,2)
|
|
virialsym(3) = virial_tmp(3,3)
|
|
virialsym(4) = virial_tmp(1,2)
|
|
virialsym(5) = virial_tmp(1,3)
|
|
virialsym(6) = virial_tmp(2,3)
|
|
do k1 = 1,6
|
|
virial(k1) = virial(k1) + virialsym(k1)
|
|
end do
|
|
|
|
if (Latomvirial.eq.1) then
|
|
frac = 1.0d0/idbo1(ncubo)
|
|
do k1 = 1,6
|
|
vtmp = virialsym(k1)*frac
|
|
do i4=1,idbo1(ncubo)
|
|
ihu=idbo(ncubo,i4)
|
|
atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
endif
|
|
|
|
end do
|
|
end do
|
|
|
|
|
|
end do
|
|
|
|
|
|
return
|
|
end
|
|
**********************************************************************
|
|
**********************************************************************
|
|
|
|
subroutine valang
|
|
|
|
**********************************************************************
|
|
#include "cbka.blk"
|
|
#include "cbkabo.blk"
|
|
#include "cbkbo.blk"
|
|
#include "cbkbopi.blk"
|
|
#include "cbkbopi2.blk"
|
|
#include "cbkconst.blk"
|
|
#include "cbkc.blk"
|
|
#include "cbkd.blk"
|
|
#include "cbkdbopi2ndc.blk"
|
|
#include "cbkdbopindc.blk"
|
|
#include "cbkdcell.blk"
|
|
#include "cbkdhdc.blk"
|
|
#include "cbkenergies.blk"
|
|
#include "cbkff.blk"
|
|
#include "cbkh.blk"
|
|
#include "cbkia.blk"
|
|
#include "cbkidbo.blk"
|
|
#include "cbklonpar.blk"
|
|
#include "cbknubon2.blk"
|
|
#include "cbkvalence.blk"
|
|
#include "control.blk"
|
|
#include "valang.blk"
|
|
#include "small.blk"
|
|
dimension j(3)
|
|
dimension virial_tmp(3,3),virialsym(6)
|
|
**********************************************************************
|
|
* *
|
|
* Calculate valency angle energies and first derivatives *
|
|
* *
|
|
**********************************************************************
|
|
c$$$ if (ndebug.eq.1) then
|
|
c$$$C open (65,file='fort.65',status='unknown',access='append')
|
|
c$$$ write (65,*) 'In valang'
|
|
c$$$ call timer(65)
|
|
c$$$ close (65)
|
|
c$$$ end if
|
|
* eco=0.0
|
|
ev=0.0
|
|
ecoa=0.0
|
|
if (nval.eq.0) return
|
|
|
|
do 10 i1=1,nval
|
|
ity=iv(i1,1)
|
|
j(1)=iv(i1,2)
|
|
j(2)=iv(i1,3)
|
|
j(3)=iv(i1,4)
|
|
|
|
if (j(2) .le. na_local) then
|
|
|
|
la=iv(i1,5)
|
|
lb=iv(i1,6)
|
|
boa=bo(la)-cutof2
|
|
bob=bo(lb)-cutof2
|
|
if (boa.lt.zero.or.bob.lt.zero) goto 10
|
|
|
|
hl=h(i1) ! Calculated earlier in routine calval
|
|
**********************************************************************
|
|
* *
|
|
* Calculate valency angle energy *
|
|
* *
|
|
**********************************************************************
|
|
nbocen=ia(j(2),2)
|
|
sbo2=0.0
|
|
vmbo=1.0
|
|
|
|
do i2=1,nbocen
|
|
ibv=nubon2(j(2),i2)
|
|
if (bo(ibv).gt.0.0) then
|
|
vmbo=vmbo*exp(-bo(ibv)**8)
|
|
sbo2=sbo2+bopi(ibv)+bopi2(ibv)
|
|
endif
|
|
end do
|
|
|
|
ity2=ia(j(2),1)
|
|
* exbo=abo(j(2))-stlp(ia(j(2),1))
|
|
exbo=abo(j(2))-valf(ity2)
|
|
* if (exbo.gt.zero) exbo=zero
|
|
* expov=exp(vka8(ity)*exbo)
|
|
* expov2=exp(-vpar(13)*exbo)
|
|
* htov1=2.0+expov2
|
|
* htov2=1.0+expov+expov2
|
|
* evboadj=htov1/htov2
|
|
evboadj=1.0
|
|
expun=exp(-vkac(ity)*exbo)
|
|
expun2=exp(vpar(15)*exbo)
|
|
htun1=2.0+expun2
|
|
htun2=1.0+expun+expun2
|
|
evboadj2=vval4(ity2)-(vval4(ity2)-1.0)*htun1/htun2
|
|
**********************************************************************
|
|
* *
|
|
* Calculate number of lone pairs *
|
|
* *
|
|
**********************************************************************
|
|
dsbo2dvlp=(1.0-vmbo)
|
|
vlpadj=zero
|
|
exlp1=abo(j(2))-stlp(ia(j(2),1))
|
|
exlp2=2.0*int(exlp1/2.0)
|
|
exlp=exlp1-exlp2
|
|
if (exlp.lt.zero) then
|
|
* expvlp=exp(-vpar(16)*(2.0+exlp)*(2.0+exlp))
|
|
* vlpadj=expvlp-int(exlp1/2.0)
|
|
* dsbo2dvlp=(1.0-vmbo)*(1.0-vpar(34)*2.0*
|
|
* $(2.0+exlp)*vpar(16)*expvlp)
|
|
vlpadj=vlp(j(2))
|
|
dsbo2dvlp=(1.0-vmbo)*(1.0+vpar(34)*dvlpdsbo(j(2)))
|
|
end if
|
|
|
|
sbo2=sbo2+(1.0-vmbo)*(-exbo-vpar(34)*vlpadj)
|
|
dsbo2dvmbo=exbo+vpar(34)*vlpadj
|
|
|
|
sbo2h=sbo2
|
|
powv=vpar(17)
|
|
if (sbo2.le.0.0) sbo2h=0.0
|
|
if (sbo2.gt.0.0.and.sbo2.le.1.0) sbo2h=sbo2**powv
|
|
if (sbo2.gt.1.0.and.sbo2.lt.2.0) sbo2h=2.0-(2.0-sbo2)**powv
|
|
if (sbo2.gt.2.0) sbo2h=2.0
|
|
thba=th0(ity)
|
|
expsbo=exp(-vpar(18)*(2.0-sbo2h))
|
|
thetao=180.0-thba*(1.0-expsbo)
|
|
|
|
thetao=thetao*dgrrdn
|
|
thdif=(thetao-hl)
|
|
thdi2=thdif*thdif
|
|
dthsbo=dgrrdn*thba*vpar(18)*expsbo
|
|
if (sbo2.lt.0.0) dthsbo=zero
|
|
if (sbo2.gt.0.0.and.sbo2.le.1.0)
|
|
$dthsbo=powv*(sbo2**(powv-1.0))*dgrrdn*thba*vpar(18)*expsbo
|
|
if (sbo2.gt.1.0.and.sbo2.lt.2.0)
|
|
$dthsbo=powv*((2.0-sbo2)**(powv-1.0))*dgrrdn*thba*vpar(18)*expsbo
|
|
if (sbo2.gt.2.0) dthsbo=zero
|
|
|
|
exphu=vka(ity)*exp(-vka3(ity)*thdi2)
|
|
exphu2=vka(ity)-exphu
|
|
if (vka(ity).lt.zero) exphu2=exphu2-vka(ity) !To avoid linear Me-H-Me angles (6/6/06)
|
|
boap=boa**vval2(ity)
|
|
boap2=boa**(vval2(ity)-1.0)
|
|
bobp=bob**vval2(ity)
|
|
bobp2=bob**(vval2(ity)-1.0)
|
|
exa=exp(-vval1(ity2)*boap)
|
|
exb=exp(-vval1(ity2)*bobp)
|
|
dexadboa=vval2(ity)*vval1(ity2)*exa*boap2
|
|
dexbdbob=vval2(ity)*vval1(ity2)*exb*bobp2
|
|
exa2=(1.0-exa)
|
|
exb2=(1.0-exb)
|
|
|
|
evh=evboadj2*evboadj*exa2*exb2*exphu2
|
|
devdlb=evboadj2*evboadj*dexbdbob*exa2*exphu2
|
|
devdla=evboadj2*evboadj*dexadboa*exb2*exphu2
|
|
devdsbo=2.0*evboadj2*evboadj*dthsbo*exa2*exb2*
|
|
$vka3(ity)*thdif*exphu
|
|
devdh=-2.0*evboadj2*evboadj*exa2*exb2*vka3(ity)*thdif*exphu
|
|
|
|
devdsbo2=
|
|
$evboadj*exa2*exb2*exphu2*(vval4(ity2)-1.0)*(-vpar(15)*expun2/htun2
|
|
$+htun1*(vpar(15)*expun2-vkac(ity)*expun)/(htun2*htun2))
|
|
|
|
* devdsbo2=-evboadj2*exa2*exb2*exphu2*(vpar(13)*expov2/htov2+
|
|
* $htov1*(vka8(ity)*expov-vpar(13)*expov2)/(htov2*htov2))+
|
|
* $evboadj*exa2*exb2*exphu2*(vpar(14)-1.0)*(-vpar(15)*expun2/htun2
|
|
* $+htun1*(vpar(15)*expun2-vkac(ity)*expun)/(htun2*htun2))
|
|
|
|
if (j(2) .le. na_local) then
|
|
ev=ev+evh
|
|
estrain(j(2))=estrain(j(2))+evh !central atom energy
|
|
endif
|
|
|
|
* write (64,'(4i8,18f8.2)')mdstep,j(1),j(2),j(3),sbo2,sbo2h,
|
|
* $thetao*rdndgr,hl*rdndgr,bo(la),bo(lb),bopi(la),
|
|
* $vlp(j(2)),exbo,vlpadj,vmbo,evh,ev,vka(ity)
|
|
**********************************************************************
|
|
* *
|
|
* Calculate penalty for two double bonds in valency angle *
|
|
* *
|
|
**********************************************************************
|
|
exbo=abo(j(2))-aval(ia(j(2),1))
|
|
expov=exp(vpar(22)*exbo)
|
|
expov2=exp(-vpar(21)*exbo)
|
|
htov1=2.0+expov2
|
|
htov2=1.0+expov+expov2
|
|
ecsboadj=htov1/htov2
|
|
exphu1=exp(-vpar(20)*(boa-2.0)*(boa-2.0))
|
|
exphu2=exp(-vpar(20)*(bob-2.0)*(bob-2.0))
|
|
|
|
epenh=vkap(ity)*ecsboadj*exphu1*exphu2
|
|
estrain(j(2))=estrain(j(2))+epenh
|
|
epen=epen+epenh
|
|
decoadboa=-2.0*vpar(20)*epenh*(boa-2.0)
|
|
decoadbob=-2.0*vpar(20)*epenh*(bob-2.0)
|
|
|
|
decdsbo2=-vkap(ity)*exphu1*exphu2*(vpar(21)*expov2/htov2+htov1*
|
|
$(vpar(22)*expov-vpar(21)*expov2)/(htov2*htov2))
|
|
**********************************************************************
|
|
* *
|
|
* Calculate valency angle conjugation energy *
|
|
* *
|
|
**********************************************************************
|
|
unda=abo(j(1))-boa
|
|
* ovb=abo(j(2))-valf(ia(j(2),1))
|
|
ovb=abo(j(2))-vval3(ia(j(2),1)) !Modification for Ru 7/6/2004
|
|
|
|
undc=abo(j(3))-bob
|
|
ba=(boa-1.50)*(boa-1.50)
|
|
bb=(bob-1.50)*(bob-1.50)
|
|
exphua=exp(-vpar(31)*ba)
|
|
exphub=exp(-vpar(31)*bb)
|
|
exphuua=exp(-vpar(39)*unda*unda)
|
|
exphuob=exp(vpar(3)*ovb)
|
|
exphuuc=exp(-vpar(39)*undc*undc)
|
|
hulpob=1.0/(1.0+exphuob)
|
|
ecoah=vka8(ity)*exphua*exphub*exphuua*exphuuc*hulpob
|
|
decodbola=-2.0*vka8(ity)*(boa-1.50)*vpar(31)*exphua*exphub
|
|
$*exphuua*exphuuc*hulpob+vpar(39)*vka8(ity)*exphua*exphub*
|
|
$exphuua*exphuuc*hulpob*2.0*unda
|
|
decodbolb=-2.0*vka8(ity)*(bob-1.50)*vpar(31)*exphua*exphub
|
|
$*exphuua*exphuuc*hulpob+vpar(39)*vka8(ity)*exphua*exphub*
|
|
$exphuua*exphuuc*hulpob*2.0*undc
|
|
decodboua=-2.0*unda*vka8(ity)*vpar(39)*exphua*exphub
|
|
$*exphuua*exphuuc*hulpob
|
|
decodbouc=-2.0*undc*vka8(ity)*vpar(39)*exphua*exphub
|
|
$*exphuua*exphuuc*hulpob
|
|
decodboob=-vka8(ity)*exphua*exphub*exphuua*exphuuc*hulpob*
|
|
$hulpob*vpar(3)*exphuob
|
|
* decodboob=zero
|
|
* decodboua=zero
|
|
* decodbouc=zero
|
|
|
|
ecoa=ecoa+ecoah
|
|
estrain(j(2))=estrain(j(2))+ecoah !central atom energy
|
|
|
|
**********************************************************************
|
|
* *
|
|
* Calculate derivative valency energy to cartesian coordinates *
|
|
* *
|
|
**********************************************************************
|
|
if (Lvirial.eq.1) then
|
|
do k1=1,3
|
|
do k2=1,3
|
|
virial_tmp(k1,k2) = 0.0
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
do k1=1,3
|
|
do k2=1,3
|
|
ftmp = devdh*dhdc(i1,k1,k2)
|
|
d(k1,j(k2))=d(k1,j(k2))+ftmp
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1p=1,3
|
|
virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(j(k2),k1p)
|
|
end do
|
|
endif
|
|
|
|
end do
|
|
end do
|
|
if (Lvirial.eq.1) then
|
|
virialsym(1) = virial_tmp(1,1)
|
|
virialsym(2) = virial_tmp(2,2)
|
|
virialsym(3) = virial_tmp(3,3)
|
|
virialsym(4) = virial_tmp(1,2)
|
|
virialsym(5) = virial_tmp(1,3)
|
|
virialsym(6) = virial_tmp(2,3)
|
|
do k1 = 1,6
|
|
virial(k1) = virial(k1) + virialsym(k1)
|
|
end do
|
|
|
|
if (Latomvirial.eq.1) then
|
|
frac = 1.0d0/3
|
|
do k1 = 1,6
|
|
vtmp = virialsym(k1)*frac
|
|
do k2=1,3
|
|
ihu=j(k2)
|
|
atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
endif
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1=1,3
|
|
do k2=1,3
|
|
virial_tmp(k1,k2) = 0.0
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
do i2=1,idbo1(la)
|
|
ihu=idbo(la,i2)
|
|
do k1=1,3
|
|
ftmp = (devdla+decoadboa+decodbola)*
|
|
$dbondc(la,k1,i2)
|
|
d(k1,ihu)=d(k1,ihu)+ftmp
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1p=1,3
|
|
virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p)
|
|
end do
|
|
endif
|
|
|
|
end do
|
|
end do
|
|
if (Lvirial.eq.1) then
|
|
virialsym(1) = virial_tmp(1,1)
|
|
virialsym(2) = virial_tmp(2,2)
|
|
virialsym(3) = virial_tmp(3,3)
|
|
virialsym(4) = virial_tmp(1,2)
|
|
virialsym(5) = virial_tmp(1,3)
|
|
virialsym(6) = virial_tmp(2,3)
|
|
do k1 = 1,6
|
|
virial(k1) = virial(k1) + virialsym(k1)
|
|
end do
|
|
|
|
if (Latomvirial.eq.1) then
|
|
frac = 1.0d0/idbo1(la)
|
|
do k1 = 1,6
|
|
vtmp = virialsym(k1)*frac
|
|
do i2=1,idbo1(la)
|
|
ihu=idbo(la,i2)
|
|
atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
endif
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1=1,3
|
|
do k2=1,3
|
|
virial_tmp(k1,k2) = 0.0
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
do i2=1,idbo1(lb)
|
|
ihu=idbo(lb,i2)
|
|
do k1=1,3
|
|
ftmp = (devdlb+decoadbob+decodbolb)*
|
|
$dbondc(lb,k1,i2)
|
|
d(k1,ihu)=d(k1,ihu)+ftmp
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1p=1,3
|
|
virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p)
|
|
end do
|
|
endif
|
|
|
|
end do
|
|
end do
|
|
if (Lvirial.eq.1) then
|
|
virialsym(1) = virial_tmp(1,1)
|
|
virialsym(2) = virial_tmp(2,2)
|
|
virialsym(3) = virial_tmp(3,3)
|
|
virialsym(4) = virial_tmp(1,2)
|
|
virialsym(5) = virial_tmp(1,3)
|
|
virialsym(6) = virial_tmp(2,3)
|
|
do k1 = 1,6
|
|
virial(k1) = virial(k1) + virialsym(k1)
|
|
end do
|
|
|
|
if (Latomvirial.eq.1) then
|
|
frac = 1.0d0/idbo1(lb)
|
|
do k1 = 1,6
|
|
vtmp = virialsym(k1)*frac
|
|
do i2=1,idbo1(lb)
|
|
ihu=idbo(lb,i2)
|
|
atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
endif
|
|
|
|
do i2=1,nbocen
|
|
j5=ia(j(2),2+i2)
|
|
ibv=nubon2(j(2),i2)
|
|
if (bo(ibv).gt.0.0) then
|
|
dvmbodbo=-vmbo*8.0*bo(ibv)**7
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1=1,3
|
|
do k2=1,3
|
|
virial_tmp(k1,k2) = 0.0
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
do i3=1,idbo1(ibv)
|
|
ihu=idbo(ibv,i3)
|
|
do k1=1,3
|
|
ftmp = (-dsbo2dvlp*devdsbo+devdsbo2+decdsbo2
|
|
$+dvmbodbo*dsbo2dvmbo*devdsbo)*
|
|
$dbondc(ibv,k1,i3)+devdsbo*(dbopindc(ibv,k1,i3)+
|
|
$dbopi2ndc(ibv,k1,i3))
|
|
d(k1,ihu)=d(k1,ihu)+ftmp
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1p=1,3
|
|
virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p)
|
|
end do
|
|
endif
|
|
|
|
end do
|
|
end do
|
|
if (Lvirial.eq.1) then
|
|
virialsym(1) = virial_tmp(1,1)
|
|
virialsym(2) = virial_tmp(2,2)
|
|
virialsym(3) = virial_tmp(3,3)
|
|
virialsym(4) = virial_tmp(1,2)
|
|
virialsym(5) = virial_tmp(1,3)
|
|
virialsym(6) = virial_tmp(2,3)
|
|
do k1 = 1,6
|
|
virial(k1) = virial(k1) + virialsym(k1)
|
|
end do
|
|
|
|
if (Latomvirial.eq.1) then
|
|
frac = 1.0d0/idbo1(ibv)
|
|
do k1 = 1,6
|
|
vtmp = virialsym(k1)*frac
|
|
do i3=1,idbo1(ibv)
|
|
ihu=idbo(ibv,i3)
|
|
atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
endif
|
|
|
|
endif
|
|
end do
|
|
|
|
do i2=1,ia(j(1),2)
|
|
j5=ia(j(1),2+i2)
|
|
ibv=nubon2(j(1),i2)
|
|
if (bo(ibv).gt.0.0) then
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1=1,3
|
|
do k2=1,3
|
|
virial_tmp(k1,k2) = 0.0
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
do i3=1,idbo1(ibv)
|
|
ihu=idbo(ibv,i3)
|
|
do k1=1,3
|
|
ftmp = decodboua*dbondc(ibv,k1,i3)
|
|
d(k1,ihu)=d(k1,ihu)+ftmp
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1p=1,3
|
|
virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p)
|
|
end do
|
|
endif
|
|
|
|
end do
|
|
end do
|
|
if (Lvirial.eq.1) then
|
|
virialsym(1) = virial_tmp(1,1)
|
|
virialsym(2) = virial_tmp(2,2)
|
|
virialsym(3) = virial_tmp(3,3)
|
|
virialsym(4) = virial_tmp(1,2)
|
|
virialsym(5) = virial_tmp(1,3)
|
|
virialsym(6) = virial_tmp(2,3)
|
|
do k1 = 1,6
|
|
virial(k1) = virial(k1) + virialsym(k1)
|
|
end do
|
|
|
|
if (Latomvirial.eq.1) then
|
|
frac = 1.0d0/idbo1(ibv)
|
|
do k1 = 1,6
|
|
vtmp = virialsym(k1)*frac
|
|
do i3=1,idbo1(ibv)
|
|
ihu=idbo(ibv,i3)
|
|
atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
endif
|
|
|
|
endif
|
|
end do
|
|
|
|
do i2=1,ia(j(2),2)
|
|
j5=ia(j(2),2+i2)
|
|
ibv=nubon2(j(2),i2)
|
|
if (bo(ibv).gt.0.0) then
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1=1,3
|
|
do k2=1,3
|
|
virial_tmp(k1,k2) = 0.0
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
do i3=1,idbo1(ibv)
|
|
ihu=idbo(ibv,i3)
|
|
do k1=1,3
|
|
ftmp = decodboob*dbondc(ibv,k1,i3)
|
|
d(k1,ihu)=d(k1,ihu)+ftmp
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1p=1,3
|
|
virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p)
|
|
end do
|
|
endif
|
|
|
|
end do
|
|
end do
|
|
if (Lvirial.eq.1) then
|
|
virialsym(1) = virial_tmp(1,1)
|
|
virialsym(2) = virial_tmp(2,2)
|
|
virialsym(3) = virial_tmp(3,3)
|
|
virialsym(4) = virial_tmp(1,2)
|
|
virialsym(5) = virial_tmp(1,3)
|
|
virialsym(6) = virial_tmp(2,3)
|
|
do k1 = 1,6
|
|
virial(k1) = virial(k1) + virialsym(k1)
|
|
end do
|
|
|
|
if (Latomvirial.eq.1) then
|
|
frac = 1.0d0/idbo1(ibv)
|
|
do k1 = 1,6
|
|
vtmp = virialsym(k1)*frac
|
|
do i3=1,idbo1(ibv)
|
|
ihu=idbo(ibv,i3)
|
|
atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
endif
|
|
|
|
endif
|
|
end do
|
|
|
|
do i2=1,ia(j(3),2)
|
|
j5=ia(j(3),2+i2)
|
|
ibv=nubon2(j(3),i2)
|
|
if (bo(ibv).gt.0.0) then
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1=1,3
|
|
do k2=1,3
|
|
virial_tmp(k1,k2) = 0.0
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
do i3=1,idbo1(ibv)
|
|
ihu=idbo(ibv,i3)
|
|
do k1=1,3
|
|
ftmp = decodbouc*dbondc(ibv,k1,i3)
|
|
d(k1,ihu)=d(k1,ihu)+ftmp
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1p=1,3
|
|
virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p)
|
|
end do
|
|
endif
|
|
|
|
end do
|
|
end do
|
|
if (Lvirial.eq.1) then
|
|
virialsym(1) = virial_tmp(1,1)
|
|
virialsym(2) = virial_tmp(2,2)
|
|
virialsym(3) = virial_tmp(3,3)
|
|
virialsym(4) = virial_tmp(1,2)
|
|
virialsym(5) = virial_tmp(1,3)
|
|
virialsym(6) = virial_tmp(2,3)
|
|
do k1 = 1,6
|
|
virial(k1) = virial(k1) + virialsym(k1)
|
|
end do
|
|
|
|
if (Latomvirial.eq.1) then
|
|
frac = 1.0d0/idbo1(ibv)
|
|
do k1 = 1,6
|
|
vtmp = virialsym(k1)*frac
|
|
do i3=1,idbo1(ibv)
|
|
ihu=idbo(ibv,i3)
|
|
atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
endif
|
|
|
|
endif
|
|
end do
|
|
|
|
endif
|
|
|
|
10 continue
|
|
|
|
return
|
|
end
|
|
**********************************************************************
|
|
**********************************************************************
|
|
|
|
subroutine hbond
|
|
|
|
**********************************************************************
|
|
#include "cbka.blk"
|
|
#include "cbkbo.blk"
|
|
#include "cbkconst.blk"
|
|
#include "cbkc.blk"
|
|
#include "cbkd.blk"
|
|
#include "cbkdcell.blk"
|
|
#include "cbkenergies.blk"
|
|
#include "cbkidbo.blk"
|
|
#include "cbksrthb.blk"
|
|
#include "control.blk"
|
|
#include "cbkhbond.blk"
|
|
#include "small.blk"
|
|
dimension drda(3),j(3),dvdc(3,3),dargdc(3,3)
|
|
dimension virial_tmp(3,3),virialsym(6)
|
|
**********************************************************************
|
|
* *
|
|
* Calculate hydrogen bond energies and first derivatives *
|
|
* *
|
|
**********************************************************************
|
|
c$$$ if (ndebug.eq.1) then
|
|
c$$$C open (65,file='fort.65',status='unknown',access='append')
|
|
c$$$ write (65,*) 'In hbond'
|
|
c$$$ call timer(65)
|
|
c$$$ close (65)
|
|
c$$$ end if
|
|
ehb=zero
|
|
do 10 i1=1,nhb
|
|
ityhb=ihb(i1,1)
|
|
j(1)=ihb(i1,2)
|
|
j(2)=ihb(i1,3)
|
|
j(3)=ihb(i1,4)
|
|
la=ihb(i1,5)
|
|
boa=bo(la)
|
|
call dista2(j(2),j(3),rda,dxm,dym,dzm)
|
|
drda(1)=dxm/rda
|
|
drda(2)=dym/rda
|
|
drda(3)=dzm/rda
|
|
call calvalhb(j(1),j(2),j(3),ix,iy,iz,arg,hhb(i1),dvdc,dargdc)
|
|
rhu1=rhb(ityhb)/rda
|
|
rhu2=rda/rhb(ityhb)
|
|
sinhu=sin(hhb(i1)/2.0)
|
|
sin2=sinhu*sinhu
|
|
exphu1=exp(-vhb1(ityhb)*boa)
|
|
exphu2=exp(-vhb2(ityhb)*(rhu1+rhu2-2.0))
|
|
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
|
|
|
|
**********************************************************************
|
|
* *
|
|
* Calculate first derivatives *
|
|
* *
|
|
**********************************************************************
|
|
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
|
|
do k2=1,3
|
|
virial_tmp(k1,k2) = 0.0
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
do k1=1,3
|
|
ftmp = dehbdrda*drda(k1)
|
|
d(k1,j(2))=d(k1,j(2))+ftmp
|
|
d(k1,j(3))=d(k1,j(3))-ftmp
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1p=1,3
|
|
virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+
|
|
$ ftmp*c(j(2),k1p)-ftmp*c(j(3),k1p)
|
|
end do
|
|
endif
|
|
|
|
end do
|
|
if (Lvirial.eq.1) then
|
|
virialsym(1) = virial_tmp(1,1)
|
|
virialsym(2) = virial_tmp(2,2)
|
|
virialsym(3) = virial_tmp(3,3)
|
|
virialsym(4) = virial_tmp(1,2)
|
|
virialsym(5) = virial_tmp(1,3)
|
|
virialsym(6) = virial_tmp(2,3)
|
|
do k1 = 1,6
|
|
virial(k1) = virial(k1) + virialsym(k1)
|
|
end do
|
|
|
|
if (Latomvirial.eq.1) then
|
|
frac = 1.0d0/2
|
|
do k1 = 1,6
|
|
vtmp = virialsym(k1)*frac
|
|
ihu = j(2)
|
|
atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp
|
|
ihu = j(3)
|
|
atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp
|
|
end do
|
|
endif
|
|
|
|
endif
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1=1,3
|
|
do k2=1,3
|
|
virial_tmp(k1,k2) = 0.0
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
do k1=1,3
|
|
do k2=1,3
|
|
ftmp = dehbdv*dvdc(k1,k2)
|
|
d(k1,j(k2))=d(k1,j(k2))+ftmp
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1p=1,3
|
|
virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(j(k2),k1p)
|
|
end do
|
|
endif
|
|
|
|
end do
|
|
end do
|
|
if (Lvirial.eq.1) then
|
|
virialsym(1) = virial_tmp(1,1)
|
|
virialsym(2) = virial_tmp(2,2)
|
|
virialsym(3) = virial_tmp(3,3)
|
|
virialsym(4) = virial_tmp(1,2)
|
|
virialsym(5) = virial_tmp(1,3)
|
|
virialsym(6) = virial_tmp(2,3)
|
|
do k1 = 1,6
|
|
virial(k1) = virial(k1) + virialsym(k1)
|
|
end do
|
|
|
|
if (Latomvirial.eq.1) then
|
|
frac = 1.0d0/3
|
|
do k1 = 1,6
|
|
vtmp = virialsym(k1)*frac
|
|
do k2=1,3
|
|
ihu=j(k2)
|
|
atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
endif
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1=1,3
|
|
do k2=1,3
|
|
virial_tmp(k1,k2) = 0.0
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
do i2=1,idbo1(la)
|
|
ihu=idbo(la,i2)
|
|
do k1=1,3
|
|
ftmp = dehbdbo*dbondc(la,k1,i2)
|
|
d(k1,ihu)=d(k1,ihu)+ftmp
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1p=1,3
|
|
virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p)
|
|
end do
|
|
endif
|
|
|
|
end do
|
|
end do
|
|
if (Lvirial.eq.1) then
|
|
virialsym(1) = virial_tmp(1,1)
|
|
virialsym(2) = virial_tmp(2,2)
|
|
virialsym(3) = virial_tmp(3,3)
|
|
virialsym(4) = virial_tmp(1,2)
|
|
virialsym(5) = virial_tmp(1,3)
|
|
virialsym(6) = virial_tmp(2,3)
|
|
do k1 = 1,6
|
|
virial(k1) = virial(k1) + virialsym(k1)
|
|
end do
|
|
|
|
if (Latomvirial.eq.1) then
|
|
frac = 1.0d0/idbo1(la)
|
|
do k1 = 1,6
|
|
vtmp = virialsym(k1)*frac
|
|
do i2=1,idbo1(la)
|
|
ihu=idbo(la,i2)
|
|
atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
endif
|
|
|
|
10 continue
|
|
return
|
|
end
|
|
|
|
**********************************************************************
|
|
**********************************************************************
|
|
|
|
subroutine torang
|
|
|
|
**********************************************************************
|
|
#include "cbka.blk"
|
|
#include "cbkabo.blk"
|
|
#include "cbkbo.blk"
|
|
#include "cbkbopi.blk"
|
|
#include "cbkc.blk"
|
|
#include "cbkconst.blk"
|
|
#include "cbkd.blk"
|
|
#include "cbkdbopindc.blk"
|
|
#include "cbkdcell.blk"
|
|
#include "cbkdhdc.blk"
|
|
#include "cbkdrdc.blk"
|
|
#include "cbkenergies.blk"
|
|
#include "cbkff.blk"
|
|
#include "cbkfftorang.blk"
|
|
#include "cbkh.blk"
|
|
#include "cbkia.blk"
|
|
#include "cbkidbo.blk"
|
|
#include "cbkinit.blk"
|
|
#include "cbknubon2.blk"
|
|
#include "cbkrbo.blk"
|
|
#include "cbktorang.blk"
|
|
#include "cbktorsion.blk"
|
|
#include "cbktregime.blk"
|
|
#include "cbkvalence.blk"
|
|
#include "cellcoord.blk"
|
|
#include "control.blk"
|
|
#include "small.blk"
|
|
|
|
DIMENSION A(3),DRDA(3),DADC(4),DRADC(3,4),DRBDC(3,4),
|
|
$DRCDC(3,4),DHDDC(3,4),DHEDC(3,4),DRVDC(3,4),DTDC(3,4),
|
|
$DNDC(3,4)
|
|
dimension j(4),dh1rdc(3,3),dh2rdc(3,3),dargdc(3,3)
|
|
dimension virial_tmp(3,3),virialsym(6)
|
|
**********************************************************************
|
|
* *
|
|
* Calculate torsion angle energies and first derivatives *
|
|
* *
|
|
**********************************************************************
|
|
c$$$ if (ndebug.eq.1) then
|
|
c$$$C open (65,file='fort.65',status='unknown',access='append')
|
|
c$$$ write (65,*) 'In torang'
|
|
c$$$ call timer(65)
|
|
c$$$ close (65)
|
|
c$$$ end if
|
|
do k1=1,3
|
|
do k2=1,4
|
|
dhddc(k1,k2)=0.0
|
|
dhedc(k1,k2)=0.0
|
|
dradc(k1,k2)=0.0
|
|
drbdc(k1,k2)=0.0
|
|
drcdc(k1,k2)=0.0
|
|
end do
|
|
end do
|
|
et=0.0
|
|
eth12=0.0
|
|
eco=0.0
|
|
dadc(1)=1.0
|
|
dadc(2)=0.0
|
|
dadc(3)=0.0
|
|
dadc(4)=-1.0
|
|
if (ntor.eq.0) return
|
|
|
|
do 10 i1=1,ntor
|
|
j(1)=it(i1,2)
|
|
j(2)=it(i1,3)
|
|
j(3)=it(i1,4)
|
|
j(4)=it(i1,5)
|
|
|
|
ity=it(i1,1)
|
|
la=it(i1,6)
|
|
lb=it(i1,7)
|
|
lc=it(i1,8)
|
|
call calvalres(j(1),j(2),j(3),arg1,ht1,dh1rdc,dargdc)
|
|
call calvalres(j(2),j(3),j(4),arg2,ht2,dh2rdc,dargdc)
|
|
boa=bo(la)-cutof2
|
|
bob=bo(lb)-cutof2
|
|
boc=bo(lc)-cutof2
|
|
if (boa.lt.zero.or.bob.lt.zero.or.boc.lt.zero)
|
|
$goto 10
|
|
r42=0.0
|
|
ivl1=ibsym(la)
|
|
ivl2=ibsym(lb)
|
|
ivl3=ibsym(lc)
|
|
isign1=1
|
|
isign2=1
|
|
isign3=1
|
|
rla=rbo(la)
|
|
rlb=rbo(lb)
|
|
|
|
call dista2(j(1),j(4),r4,a(1),a(2),a(3))
|
|
**********************************************************************
|
|
* *
|
|
* Determine torsion angle *
|
|
* *
|
|
**********************************************************************
|
|
d142=r4*r4
|
|
rla=rbo(la)
|
|
rlb=rbo(lb)
|
|
rlc=rbo(lc)
|
|
coshd=cos(ht1)
|
|
coshe=cos(ht2)
|
|
sinhd=sin(ht1)
|
|
sinhe=sin(ht2)
|
|
poem=2.0*rla*rlc*sinhd*sinhe
|
|
poem2=poem*poem
|
|
tel=rla*rla+rlb*rlb+rlc*rlc-d142-2.0*(rla*rlb*coshd-rla*rlc*
|
|
$coshd*coshe+rlb*rlc*coshe)
|
|
if (poem.lt.1e-20) poem=1e-20
|
|
arg=tel/poem
|
|
if (arg.gt.1.0) arg=1.0
|
|
if (arg.lt.-1.0) arg=-1.0
|
|
arg2=arg*arg
|
|
thg(i1)=acos(arg)*rdndgr
|
|
k1=j(1)
|
|
k2=j(2)
|
|
k3=j(3)
|
|
k4=j(4)
|
|
call dista2(k3,k2,dis,x3,y3,z3)
|
|
y32z32=y3*y3+z3*z3
|
|
wort1=sqrt(y32z32)+1e-6
|
|
wort2=sqrt(y32z32+x3*x3)+1e-6
|
|
* if (wort1.lt.1e-6) wort1=1e-6
|
|
* if (wort2.lt.1e-6) wort2=1e-6
|
|
sinalf=y3/wort1
|
|
cosalf=z3/wort1
|
|
sinbet=x3/wort2
|
|
cosbet=wort1/wort2
|
|
call dista2(k1,k2,dis,x1,y1,z1)
|
|
x1=x1*cosbet-y1*sinalf*sinbet-z1*cosalf*sinbet
|
|
y1=y1*cosalf-z1*sinalf
|
|
wort3=sqrt(x1*x1+y1*y1)+1e-6
|
|
* if (wort3.lt.1e-6) wort3=1e-6
|
|
singam=y1/wort3
|
|
cosgam=x1/wort3
|
|
call dista2(k4,k2,dis,x4,y4,z4)
|
|
x4=x4*cosbet-y4*sinalf*sinbet-z4*cosalf*sinbet
|
|
y4=y4*cosalf-z4*sinalf
|
|
y4=x4*singam-y4*cosgam
|
|
if (y4.gt.0.0) thg(i1)=-thg(i1)
|
|
if (thg(i1).lt.-179.999999d0) thg(i1)=-179.999999d0
|
|
if (thg(i1).gt.179.999999d0) thg(i1)=179.999999d0
|
|
th2=thg(i1)*dgrrdn
|
|
**********************************************************************
|
|
* *
|
|
* Calculate torsion angle energy *
|
|
* *
|
|
**********************************************************************
|
|
exbo1=abo(j(2))-valf(ia(j(2),1))
|
|
exbo2=abo(j(3))-valf(ia(j(3),1))
|
|
htovt=exbo1+exbo2
|
|
expov=exp(vpar(26)*htovt)
|
|
expov2=exp(-vpar(25)*(htovt))
|
|
htov1=2.0+expov2
|
|
htov2=1.0+expov+expov2
|
|
etboadj=htov1/htov2
|
|
|
|
btb2=bopi(lb)-1.0+etboadj
|
|
bo2t=1.0-btb2
|
|
bo2p=bo2t*bo2t
|
|
bocor2=exp(v4(ity)*bo2p)
|
|
|
|
hsin=sinhd*sinhe
|
|
ethhulp=0.50*v1(ity)*(1.0+arg)+v2(ity)*bocor2*(1.0-arg2)+
|
|
$v3(ity)*(0.50+2.0*arg2*arg-1.50*arg)
|
|
|
|
exphua=exp(-vpar(24)*boa)
|
|
exphub=exp(-vpar(24)*bob)
|
|
exphuc=exp(-vpar(24)*boc)
|
|
bocor4=(1.0-exphua)*(1.0-exphub)*(1.0-exphuc)
|
|
eth=hsin*ethhulp*bocor4
|
|
|
|
detdar=hsin*bocor4*(0.50*v1(ity)-2.0*v2(ity)*bocor2*arg+
|
|
$v3(ity)*(6.0*arg2-1.5d0))
|
|
detdhd=coshd*sinhe*bocor4*ethhulp
|
|
detdhe=sinhd*coshe*bocor4*ethhulp
|
|
|
|
detdboa=vpar(24)*exphua*(1.0-exphub)*(1.0-exphuc)*ethhulp*hsin
|
|
detdbopib=-bocor4*2.0*v4(ity)*v2(ity)*
|
|
$bo2t*bocor2*(1.0-arg2)*hsin
|
|
detdbob=vpar(24)*exphub*(1.0-exphua)*
|
|
$(1.0-exphuc)*ethhulp*hsin
|
|
detdboc=vpar(24)*exphuc*(1.0-exphua)*
|
|
$(1.0-exphub)*ethhulp*hsin
|
|
|
|
detdsbo1=-(detdbopib)*
|
|
$(vpar(25)*expov2/htov2+htov1*
|
|
$(vpar(26)*expov-vpar(25)*expov2)/(htov2*htov2))
|
|
|
|
et=et+eth
|
|
estrain(j(2))=estrain(j(2))+0.50*eth !2nd atom energy
|
|
estrain(j(3))=estrain(j(3))+0.50*eth !3rd atom energy
|
|
|
|
**********************************************************************
|
|
* *
|
|
* Calculate conjugation energy *
|
|
* *
|
|
**********************************************************************
|
|
ba=(boa-1.50)*(boa-1.50)
|
|
bb=(bob-1.50)*(bob-1.50)
|
|
bc=(boc-1.50)*(boc-1.50)
|
|
exphua1=exp(-vpar(28)*ba)
|
|
exphub1=exp(-vpar(28)*bb)
|
|
exphuc1=exp(-vpar(28)*bc)
|
|
sbo=exphua1*exphub1*exphuc1
|
|
dbohua=-2.0*(boa-1.50)*vpar(28)*exphua1*exphub1*exphuc1
|
|
dbohub=-2.0*(bob-1.50)*vpar(28)*exphua1*exphub1*exphuc1
|
|
dbohuc=-2.0*(boc-1.50)*vpar(28)*exphua1*exphub1*exphuc1
|
|
arghu0=(arg2-1.0)*sinhd*sinhe
|
|
ehulp=vconj(ity)*(arghu0+1.0)
|
|
ecoh=ehulp*sbo
|
|
decodar=sbo*vconj(ity)*2.0*arg*sinhd*sinhe
|
|
decodbola=dbohua*ehulp
|
|
decodbolb=dbohub*ehulp
|
|
decodbolc=dbohuc*ehulp
|
|
decodhd=coshd*sinhe*vconj(ity)*sbo*(arg2-1.0)
|
|
decodhe=coshe*sinhd*vconj(ity)*sbo*(arg2-1.0)
|
|
|
|
eco=eco+ecoh
|
|
estrain(j(2))=estrain(j(2))+0.50*ecoh !2nd atom energy
|
|
estrain(j(3))=estrain(j(3))+0.50*ecoh !3rd atom energy
|
|
|
|
1 continue
|
|
**********************************************************************
|
|
* *
|
|
* Calculate derivative torsion angle and conjugation energy *
|
|
* to cartesian coordinates *
|
|
* *
|
|
**********************************************************************
|
|
SINTH=SIN(THG(i1)*DGRRDN)
|
|
IF (SINTH.GE.0.0.AND.SINTH.LT.1.0D-10) SINTH=1.0D-10
|
|
IF (SINTH.LT.0.0.AND.SINTH.GT.-1.0D-10) SINTH=-1.0D-10
|
|
IF (j(1).EQ.IB(LA,2)) THEN
|
|
DO K1=1,3
|
|
DRADC(K1,1)=DRDC(LA,K1,1)
|
|
DRADC(K1,2)=DRDC(LA,K1,2)
|
|
end do
|
|
ELSE
|
|
DO K1=1,3
|
|
DRADC(K1,1)=DRDC(LA,K1,2)
|
|
DRADC(K1,2)=DRDC(LA,K1,1)
|
|
end do
|
|
ENDIF
|
|
IF (j(2).EQ.IB(LB,2)) THEN
|
|
DO K1=1,3
|
|
DRBDC(K1,2)=DRDC(LB,K1,1)
|
|
DRBDC(K1,3)=DRDC(LB,K1,2)
|
|
end do
|
|
ELSE
|
|
DO K1=1,3
|
|
DRBDC(K1,2)=DRDC(LB,K1,2)
|
|
DRBDC(K1,3)=DRDC(LB,K1,1)
|
|
end do
|
|
ENDIF
|
|
IF (j(3).EQ.IB(LC,2)) THEN
|
|
DO K1=1,3
|
|
DRCDC(K1,3)=DRDC(LC,K1,1)
|
|
DRCDC(K1,4)=DRDC(LC,K1,2)
|
|
end do
|
|
ELSE
|
|
DO K1=1,3
|
|
DRCDC(K1,3)=DRDC(LC,K1,2)
|
|
DRCDC(K1,4)=DRDC(LC,K1,1)
|
|
end do
|
|
ENDIF
|
|
|
|
do k1=1,3
|
|
dhddc(1,k1)=dh1rdc(1,k1)
|
|
dhddc(2,k1)=dh1rdc(2,k1)
|
|
dhddc(3,k1)=dh1rdc(3,k1)
|
|
dhedc(1,k1+1)=dh2rdc(1,k1)
|
|
dhedc(2,k1+1)=dh2rdc(2,k1)
|
|
dhedc(3,k1+1)=dh2rdc(3,k1)
|
|
end do
|
|
|
|
**********************************************************************
|
|
* write (64,*)j(1),j(2),j(3),j(4)
|
|
* do k1=1,3
|
|
* write (64,'(10f12.4)')(dh1rdc(k1,k2),k2=1,3),
|
|
* $(dhdc(ld,k1,k2),k2=1,3),(dhddc(k1,k2),k2=1,4)
|
|
* write (64,'(10f12.4)')(dh2rdc(k1,k2),k2=1,3),
|
|
* $(dhdc(le,k1,k2),k2=1,3),(dhedc(k1,k2),k2=1,4)
|
|
* end do
|
|
* write (64,*)
|
|
**********************************************************************
|
|
HTRA=RLA+COSHD*(RLC*COSHE-RLB)
|
|
HTRB=RLB-RLA*COSHD-RLC*COSHE
|
|
HTRC=RLC+COSHE*(RLA*COSHD-RLB)
|
|
HTHD=RLA*SINHD*(RLB-RLC*COSHE)
|
|
HTHE=RLC*SINHE*(RLB-RLA*COSHD)
|
|
HNRA=RLC*SINHD*SINHE
|
|
HNRC=RLA*SINHD*SINHE
|
|
HNHD=RLA*RLC*COSHD*SINHE
|
|
HNHE=RLA*RLC*SINHD*COSHE
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1=1,3
|
|
do k2=1,3
|
|
virial_tmp(k1,k2) = 0.0
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
DO K1=1,3
|
|
DRDA(K1)=A(K1)/R4
|
|
DO K2=1,4
|
|
DRVDC(K1,K2)=DRDA(K1)*DADC(K2)
|
|
DTDC(K1,K2)=2.0*(DRADC(K1,K2)*HTRA+DRBDC(K1,K2)*HTRB+DRCDC(K1,K2
|
|
$)*HTRC-DRVDC(K1,K2)*R4+DHDDC(K1,K2)*HTHD+DHEDC(K1,K2)*HTHE)
|
|
DNDC(K1,K2)=2.0*(DRADC(K1,K2)*HNRA+DRCDC(K1,K2)*HNRC+DHDDC(K1,K2
|
|
$)*HNHD+DHEDC(K1,K2)*HNHE)
|
|
DARGTDC(i1,K1,K2)=(DTDC(K1,K2)-ARG*DNDC(K1,K2))/POEM
|
|
|
|
ftmp = DARGTDC(i1,K1,K2)*detdar+
|
|
$dargtdc(i1,k1,k2)*decodar+(detdhd+decodhd)*dhddc(k1,k2)+
|
|
$(detdhe+decodhe)*dhedc(k1,k2)
|
|
D(K1,J(K2))=D(K1,J(K2))+ftmp
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1p=1,3
|
|
virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(j(k2),k1p)
|
|
end do
|
|
endif
|
|
|
|
end do
|
|
end do
|
|
if (Lvirial.eq.1) then
|
|
virialsym(1) = virial_tmp(1,1)
|
|
virialsym(2) = virial_tmp(2,2)
|
|
virialsym(3) = virial_tmp(3,3)
|
|
virialsym(4) = virial_tmp(1,2)
|
|
virialsym(5) = virial_tmp(1,3)
|
|
virialsym(6) = virial_tmp(2,3)
|
|
do k1 = 1,6
|
|
virial(k1) = virial(k1) + virialsym(k1)
|
|
end do
|
|
|
|
if (Latomvirial.eq.1) then
|
|
frac = 1.0d0/4
|
|
do k1 = 1,6
|
|
vtmp = virialsym(k1)*frac
|
|
do k2=1,4
|
|
ihu=j(k2)
|
|
atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
endif
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1=1,3
|
|
do k2=1,3
|
|
virial_tmp(k1,k2) = 0.0
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
do i2=1,idbo1(la)
|
|
ihu=idbo(la,i2)
|
|
do k1=1,3
|
|
ftmp = dbondc(la,k1,i2)*(detdboa+decodbola)
|
|
d(k1,ihu)=d(k1,ihu)+ftmp
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1p=1,3
|
|
virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p)
|
|
end do
|
|
endif
|
|
|
|
end do
|
|
end do
|
|
if (Lvirial.eq.1) then
|
|
virialsym(1) = virial_tmp(1,1)
|
|
virialsym(2) = virial_tmp(2,2)
|
|
virialsym(3) = virial_tmp(3,3)
|
|
virialsym(4) = virial_tmp(1,2)
|
|
virialsym(5) = virial_tmp(1,3)
|
|
virialsym(6) = virial_tmp(2,3)
|
|
do k1 = 1,6
|
|
virial(k1) = virial(k1) + virialsym(k1)
|
|
end do
|
|
|
|
if (Latomvirial.eq.1) then
|
|
frac = 1.0d0/idbo1(la)
|
|
do k1 = 1,6
|
|
vtmp = virialsym(k1)*frac
|
|
do i2=1,idbo1(la)
|
|
ihu=idbo(la,i2)
|
|
atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
endif
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1=1,3
|
|
do k2=1,3
|
|
virial_tmp(k1,k2) = 0.0
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
do i2=1,idbo1(lb)
|
|
ihu=idbo(lb,i2)
|
|
do k1=1,3
|
|
ftmp = dbondc(lb,k1,i2)*(detdbob+decodbolb)
|
|
$ +dbopindc(lb,k1,i2)*detdbopib
|
|
d(k1,ihu)=d(k1,ihu)+ftmp
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1p=1,3
|
|
virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p)
|
|
end do
|
|
endif
|
|
|
|
end do
|
|
end do
|
|
if (Lvirial.eq.1) then
|
|
virialsym(1) = virial_tmp(1,1)
|
|
virialsym(2) = virial_tmp(2,2)
|
|
virialsym(3) = virial_tmp(3,3)
|
|
virialsym(4) = virial_tmp(1,2)
|
|
virialsym(5) = virial_tmp(1,3)
|
|
virialsym(6) = virial_tmp(2,3)
|
|
do k1 = 1,6
|
|
virial(k1) = virial(k1) + virialsym(k1)
|
|
end do
|
|
|
|
if (Latomvirial.eq.1) then
|
|
frac = 1.0d0/idbo1(lb)
|
|
do k1 = 1,6
|
|
vtmp = virialsym(k1)*frac
|
|
do i2=1,idbo1(lb)
|
|
ihu=idbo(lb,i2)
|
|
atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
endif
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1=1,3
|
|
do k2=1,3
|
|
virial_tmp(k1,k2) = 0.0
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
do i2=1,idbo1(lc)
|
|
ihu=idbo(lc,i2)
|
|
do k1=1,3
|
|
ftmp = dbondc(lc,k1,i2)*(detdboc+decodbolc)
|
|
d(k1,ihu)=d(k1,ihu)+ftmp
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1p=1,3
|
|
virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p)
|
|
end do
|
|
endif
|
|
|
|
end do
|
|
end do
|
|
if (Lvirial.eq.1) then
|
|
virialsym(1) = virial_tmp(1,1)
|
|
virialsym(2) = virial_tmp(2,2)
|
|
virialsym(3) = virial_tmp(3,3)
|
|
virialsym(4) = virial_tmp(1,2)
|
|
virialsym(5) = virial_tmp(1,3)
|
|
virialsym(6) = virial_tmp(2,3)
|
|
do k1 = 1,6
|
|
virial(k1) = virial(k1) + virialsym(k1)
|
|
end do
|
|
|
|
if (Latomvirial.eq.1) then
|
|
frac = 1.0d0/idbo1(lc)
|
|
do k1 = 1,6
|
|
vtmp = virialsym(k1)*frac
|
|
do i2=1,idbo1(lc)
|
|
ihu=idbo(lc,i2)
|
|
atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
endif
|
|
|
|
do i2=1,ia(j(2),2)
|
|
iob=ia(j(2),2+i2)
|
|
ncubo=nubon2(j(2),i2)
|
|
if (bo(ncubo).gt.0.0) then
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1=1,3
|
|
do k2=1,3
|
|
virial_tmp(k1,k2) = 0.0
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
do i3=1,idbo1(ncubo)
|
|
ihu=idbo(ncubo,i3)
|
|
do k1=1,3
|
|
ftmp = detdsbo1*dbondc(ncubo,k1,i3)
|
|
d(k1,ihu)=d(k1,ihu)+ftmp
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1p=1,3
|
|
virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p)
|
|
end do
|
|
endif
|
|
|
|
end do
|
|
end do
|
|
if (Lvirial.eq.1) then
|
|
virialsym(1) = virial_tmp(1,1)
|
|
virialsym(2) = virial_tmp(2,2)
|
|
virialsym(3) = virial_tmp(3,3)
|
|
virialsym(4) = virial_tmp(1,2)
|
|
virialsym(5) = virial_tmp(1,3)
|
|
virialsym(6) = virial_tmp(2,3)
|
|
do k1 = 1,6
|
|
virial(k1) = virial(k1) + virialsym(k1)
|
|
end do
|
|
|
|
if (Latomvirial.eq.1) then
|
|
frac = 1.0d0/idbo1(ncubo)
|
|
do k1 = 1,6
|
|
vtmp = virialsym(k1)*frac
|
|
do i3=1,idbo1(ncubo)
|
|
ihu=idbo(ncubo,i3)
|
|
atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
endif
|
|
|
|
endif
|
|
end do
|
|
|
|
do i2=1,ia(j(3),2)
|
|
iob=ia(j(3),2+i2)
|
|
ncubo=nubon2(j(3),i2)
|
|
if (bo(ncubo).gt.0.0) then
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1=1,3
|
|
do k2=1,3
|
|
virial_tmp(k1,k2) = 0.0
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
do i3=1,idbo1(ncubo)
|
|
ihu=idbo(ncubo,i3)
|
|
do k1=1,3
|
|
ftmp = detdsbo1*dbondc(ncubo,k1,i3)
|
|
d(k1,ihu)=d(k1,ihu)+ftmp
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1p=1,3
|
|
virial_tmp(k1,k1p)=virial_tmp(k1,k1p)+ftmp*c(ihu,k1p)
|
|
end do
|
|
endif
|
|
|
|
end do
|
|
end do
|
|
|
|
if (Lvirial.eq.1) then
|
|
virialsym(1) = virial_tmp(1,1)
|
|
virialsym(2) = virial_tmp(2,2)
|
|
virialsym(3) = virial_tmp(3,3)
|
|
virialsym(4) = virial_tmp(1,2)
|
|
virialsym(5) = virial_tmp(1,3)
|
|
virialsym(6) = virial_tmp(2,3)
|
|
do k1 = 1,6
|
|
virial(k1) = virial(k1) + virialsym(k1)
|
|
end do
|
|
|
|
if (Latomvirial.eq.1) then
|
|
frac = 1.0d0/idbo1(ncubo)
|
|
do k1 = 1,6
|
|
vtmp = virialsym(k1)*frac
|
|
do i3=1,idbo1(ncubo)
|
|
ihu=idbo(ncubo,i3)
|
|
atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
endif
|
|
|
|
endif
|
|
end do
|
|
|
|
10 continue
|
|
|
|
return
|
|
end
|
|
**********************************************************************
|
|
**********************************************************************
|
|
|
|
subroutine nonbon
|
|
|
|
**********************************************************************
|
|
#include "cbka.blk"
|
|
#include "cbkc.blk"
|
|
#include "cbkch.blk"
|
|
#include "cbkconst.blk"
|
|
#include "cbkd.blk"
|
|
#include "cbkdcell.blk"
|
|
#include "cbkenergies.blk"
|
|
#include "cbkff.blk"
|
|
#include "cbkia.blk"
|
|
#include "cbknonbon.blk"
|
|
#include "cbkpairs.blk"
|
|
#include "cbknvlown.blk"
|
|
#include "cellcoord.blk"
|
|
#include "control.blk"
|
|
#include "small.blk"
|
|
|
|
dimension a(3),da(6)
|
|
dimension virial_tmp(3,3),virialsym(6)
|
|
**********************************************************************
|
|
* *
|
|
* Calculate vdWaals and Coulomb energies and derivatives *
|
|
* *
|
|
**********************************************************************
|
|
c$$$ if (ndebug.eq.1) then
|
|
c$$$C open (65,file='fort.65',status='unknown',access='append')
|
|
c$$$ write (65,*) 'In nonbon'
|
|
c$$$ call timer(65)
|
|
c$$$ end if
|
|
|
|
ew=0.0
|
|
ep=0.0
|
|
|
|
c1c=332.0638
|
|
third=one/three
|
|
fothird=4.0/3.0
|
|
twothird=2.0/3.0
|
|
h15=(vpar(29)-1.0)/vpar(29)
|
|
|
|
nptmp=0
|
|
nstmp=0
|
|
do 10 ivl=1,nvpair-nvlself
|
|
c Use precomputed midpoint criterion to decide if interaction is owned.
|
|
if (nvlown(ivl).eq.1) then
|
|
|
|
i1=nvl1(ivl)
|
|
i2=nvl2(ivl)
|
|
|
|
call dista2(i1,i2,rr,a(1),a(2),a(3))
|
|
if (rr.gt.swb.or.rr.lt.0.001) goto 10
|
|
|
|
ity1=ia(i1,1)
|
|
ity2=ia(i2,1)
|
|
imol1=iag(i1,3+mbond)
|
|
imol2=iag(i2,3+mbond)
|
|
rr2=rr*rr
|
|
|
|
sw=1.0
|
|
sw1=0.0
|
|
call taper(rr,rr2)
|
|
**********************************************************************
|
|
* *
|
|
* Calculate vdWaals energy *
|
|
* *
|
|
**********************************************************************
|
|
p1=p1co(ity1,ity2)
|
|
p2=p2co(ity1,ity2)
|
|
p3=p3co(ity1,ity2)
|
|
hulpw=(rr**vpar(29)+gamwco(ity1,ity2))
|
|
rrw=hulpw**(1.0/vpar(29))
|
|
h1=exp(p3*(1.0-rrw/p1))
|
|
h2=exp(0.50*p3*(1.0-rrw/p1))
|
|
|
|
ewh=p2*(h1-2.0*h2)
|
|
rrhuw=rr**(vpar(29)-1.0)
|
|
dewdr=(p2*p3/p1)*(h2-h1)*rrhuw*(hulpw**(-h15))
|
|
|
|
**********************************************************************
|
|
* *
|
|
* Calculate Coulomb energy *
|
|
* *
|
|
**********************************************************************
|
|
q1q2=ch(i1)*ch(i2)
|
|
hulp1=(rr2*rr+gamcco(ity1,ity2))
|
|
eph=c1c*q1q2/(hulp1**third)
|
|
depdr=-c1c*q1q2*rr2/(hulp1**fothird)
|
|
**********************************************************************
|
|
* *
|
|
* Taper correction *
|
|
* *
|
|
**********************************************************************
|
|
ephtap=eph*sw
|
|
depdrtap=depdr*sw+eph*sw1
|
|
ewhtap=ewh*sw
|
|
dewdrtap=dewdr*sw+ewh*sw1
|
|
|
|
* write (64,*)i1,i2,p1,p2,p3,gamwco(ity1,ity2),vpar(29),ewh,ew
|
|
ew=ew+ewhtap
|
|
ep=ep+ephtap
|
|
estrain(i1)=estrain(i1)+0.50*(ewhtap+ephtap) !1st atom energy
|
|
estrain(i2)=estrain(i2)+0.50*(ewhtap+ephtap) !2nd atom energy
|
|
|
|
**********************************************************************
|
|
* *
|
|
* Calculate derivatives vdWaals energy to cartesian *
|
|
* coordinates *
|
|
* *
|
|
**********************************************************************
|
|
|
|
if (Lvirial.eq.1) then
|
|
do k1=1,3
|
|
do k2=1,3
|
|
virial_tmp(k1,k2) = 0.0
|
|
end do
|
|
end do
|
|
endif
|
|
|
|
do k4=1,3
|
|
ftmp = (dewdrtap+depdrtap)*(a(k4)/rr)
|
|
d(k4,i1)=d(k4,i1)+ftmp
|
|
d(k4,i2)=d(k4,i2)-ftmp
|
|
if (Lvirial.eq.1) then
|
|
do k1p=1,3
|
|
virial_tmp(k4,k1p)=virial_tmp(k4,k1p)+
|
|
$ ftmp*c(i1,k1p)-ftmp*c(i2,k1p)
|
|
end do
|
|
endif
|
|
end do
|
|
|
|
if (Lvirial.eq.1) then
|
|
virialsym(1) = virial_tmp(1,1)
|
|
virialsym(2) = virial_tmp(2,2)
|
|
virialsym(3) = virial_tmp(3,3)
|
|
virialsym(4) = virial_tmp(1,2)
|
|
virialsym(5) = virial_tmp(1,3)
|
|
virialsym(6) = virial_tmp(2,3)
|
|
do k1 = 1,6
|
|
virial(k1) = virial(k1) + virialsym(k1)
|
|
end do
|
|
|
|
if (Latomvirial.eq.1) then
|
|
frac = 1.0d0/2
|
|
do k1 = 1,6
|
|
vtmp = virialsym(k1)*frac
|
|
ihu=i1
|
|
atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp
|
|
ihu=i2
|
|
atomvirial(k1,ihu) = atomvirial(k1,ihu) + vtmp
|
|
end do
|
|
endif
|
|
|
|
endif
|
|
|
|
endif
|
|
|
|
10 continue
|
|
|
|
return
|
|
end
|
|
|
|
**********************************************************************
|
|
**********************************************************************
|
|
|
|
subroutine efield
|
|
|
|
**********************************************************************
|
|
#include "cbka.blk"
|
|
#include "cbkc.blk"
|
|
#include "cbkch.blk"
|
|
#include "cbkcha.blk"
|
|
#include "cbkd.blk"
|
|
#include "cbkefield.blk"
|
|
#include "cbkenergies.blk"
|
|
#include "cbktregime.blk"
|
|
#include "control.blk"
|
|
#include "small.blk"
|
|
c$$$ if (ndebug.eq.1) then
|
|
c$$$C open (65,file='fort.65',status='unknown',access='append')
|
|
c$$$ write (65,*) 'In efield'
|
|
c$$$ call timer(65)
|
|
c$$$ close (65)
|
|
c$$$ end if
|
|
**********************************************************************
|
|
* *
|
|
* Electric field *
|
|
* *
|
|
**********************************************************************
|
|
efi=0.0
|
|
efix=0.0
|
|
efiy=0.0
|
|
efiz=0.0
|
|
c1c=332.0638 !Coulomb energy conversion
|
|
c1=23.02 !conversion from kcal to eV
|
|
|
|
if (ifieldx.eq.1) then
|
|
do i1=1,na
|
|
efih=vfieldx*c1*c1c*ch(i1)*c(i1,1)
|
|
efix=efix+efih
|
|
estrain(i1)=estrain(i1)+efih !atom energy
|
|
|
|
defidc=c1*c1c*vfieldx*ch(i1)
|
|
d(1,i1)=d(1,i1)+defidc
|
|
end do
|
|
end if
|
|
|
|
if (ifieldy.eq.1) then
|
|
do i1=1,na
|
|
efih=vfieldy*c1*c1c*ch(i1)*c(i1,2)
|
|
efiy=efiy+efih
|
|
estrain(i1)=estrain(i1)+efih !atom energy
|
|
|
|
defidc=c1*c1c*vfieldy*ch(i1)
|
|
d(2,i1)=d(2,i1)+defidc
|
|
end do
|
|
end if
|
|
|
|
if (ifieldz.eq.1) then
|
|
do i1=1,na
|
|
efih=vfieldz*c1*c1c*ch(i1)*c(i1,3)
|
|
efiz=efiz+efih
|
|
estrain(i1)=estrain(i1)+efih !atom energy
|
|
|
|
defidc=c1*c1c*vfieldz*ch(i1)
|
|
d(3,i1)=d(3,i1)+defidc
|
|
end do
|
|
end if
|
|
|
|
efi=efix+efiy+efiz
|
|
return
|
|
end
|
|
**********************************************************************
|
|
**********************************************************************
|
|
|
|
subroutine restraint
|
|
|
|
**********************************************************************
|
|
#include "cbka.blk"
|
|
#include "cbkatomcoord.blk"
|
|
#include "cbkc.blk"
|
|
#include "cbkconst.blk"
|
|
#include "cbkd.blk"
|
|
#include "cbkenergies.blk"
|
|
#include "cbkrestr.blk"
|
|
#include "cbktorang.blk"
|
|
#include "cbktorsion.blk"
|
|
#include "cbktregime.blk"
|
|
#include "control.blk"
|
|
#include "small.blk"
|
|
#include "cbkinit.blk"
|
|
dimension drda(3),j(4),dhrdc(3,3),dargdc(3,3)
|
|
**********************************************************************
|
|
* *
|
|
* Calculate distance restraint energy *
|
|
* *
|
|
**********************************************************************
|
|
c$$$ if (ndebug.eq.1) then
|
|
c$$$C open (65,file='fort.65',status='unknown',access='append')
|
|
c$$$ write (65,*) 'In restraint'
|
|
c$$$ call timer(65)
|
|
c$$$ close (65)
|
|
c$$$ end if
|
|
do i1=1,nrestra
|
|
ih1=irstra(i1,1)
|
|
ih2=irstra(i1,2)
|
|
if (itend(i1).eq.0.or.(mdstep.gt.itstart(i1).and.mdstep.lt.
|
|
$itend(i1))) then
|
|
call dista2(ih1,ih2,rr,dx,dy,dz)
|
|
diffr=rr-rrstra(i1)
|
|
* diffr=rrstra(i1)
|
|
exphu=exp(-vkrst2(i1)*(diffr*diffr))
|
|
erh=vkrstr(i1)*(1.0-exphu)
|
|
deresdr=2.0*vkrst2(i1)*diffr*vkrstr(i1)*exphu
|
|
* deresdr=-2.0*vkrst2(i1)*diffr*vkrstr(i1)*exphu
|
|
eres=eres+erh
|
|
drda(1)=dx/rr
|
|
drda(2)=dy/rr
|
|
drda(3)=dz/rr
|
|
do k1=1,3
|
|
d(k1,ih1)=d(k1,ih1)+deresdr*drda(k1)
|
|
d(k1,ih2)=d(k1,ih2)-deresdr*drda(k1)
|
|
end do
|
|
end if
|
|
end do
|
|
|
|
**********************************************************************
|
|
* *
|
|
* Calculate angle restraint energy *
|
|
* *
|
|
**********************************************************************
|
|
do i1=1,nrestrav
|
|
j(1)=irstrav(i1,1)
|
|
j(2)=irstrav(i1,2)
|
|
j(3)=irstrav(i1,3)
|
|
ittr=0
|
|
* do i2=1,nval
|
|
* if (j(1).eq.iv(i2,2).and.j(2).eq.iv(i2,3).and.j(3).eq.iv(i2,4))
|
|
* $ittr=i2
|
|
* end do
|
|
* if (ittr.eq.0) stop 'Wrong valence angle restraint'
|
|
call calvalres(j(1),j(2),j(3),arg,hr,dhrdc,dargdc)
|
|
vaval=hr*rdndgr
|
|
diffv=-(vaval-vrstra(i1))*dgrrdn
|
|
exphu=exp(-vkr2v(i1)*(diffv*diffv))
|
|
erh=vkrv(i1)*(1.0-exphu)
|
|
deresdv=-2.0*vkr2v(i1)*diffv*vkrv(i1)*exphu
|
|
eres=eres+erh
|
|
do k1=1,3
|
|
do k2=1,3
|
|
d(k1,j(k2))=d(k1,j(k2))+deresdv*dhrdc(k1,k2)
|
|
end do
|
|
end do
|
|
|
|
end do
|
|
|
|
**********************************************************************
|
|
* *
|
|
* Calculate torsion restraint energy *
|
|
* *
|
|
**********************************************************************
|
|
do i1=1,nrestrat
|
|
j(1)=irstrat(i1,1)
|
|
j(2)=irstrat(i1,2)
|
|
j(3)=irstrat(i1,3)
|
|
j(4)=irstrat(i1,4)
|
|
ittr=0
|
|
do i2=1,ntor
|
|
if (j(1).eq.it(i2,2).and.j(2).eq.it(i2,3).and.j(3).eq.it(i2,4)
|
|
$.and.j(4).eq.it(i2,5)) ittr=i2
|
|
if (j(4).eq.it(i2,2).and.j(3).eq.it(i2,3).and.j(2).eq.it(i2,4)
|
|
$.and.j(1).eq.it(i2,5)) ittr=i2
|
|
end do
|
|
if (ittr.eq.0) then
|
|
write (*,*)'Wrong torsion restraint'
|
|
write (*,*)i1,j(1),j(2),j(3),j(4)
|
|
stop 'Wrong torsion restraint'
|
|
end if
|
|
vtor=thg(ittr)
|
|
difft=-(vtor-trstra(i1))*dgrrdn
|
|
exphu=exp(-vkr2t(i1)*(difft*difft))
|
|
erh=vkrt(i1)*(1.0-exphu)
|
|
deresdt=2.0*vkr2t(i1)*difft*vkrt(i1)*exphu
|
|
if (vtor.lt.zero) deresdt=-deresdt
|
|
eres=eres+erh
|
|
do k1=1,3
|
|
do k2=1,4
|
|
d(k1,j(k2))=d(k1,j(k2))+deresdt*dargtdc(ittr,k1,k2)
|
|
end do
|
|
end do
|
|
|
|
end do
|
|
**********************************************************************
|
|
* *
|
|
* Calculate mass centre restraint energy *
|
|
* *
|
|
**********************************************************************
|
|
do i1=1,nrestram
|
|
j1=irstram(i1,2)
|
|
j2=irstram(i1,3)
|
|
j3=irstram(i1,4)
|
|
j4=irstram(i1,5)
|
|
kdir=irstram(i1,1)
|
|
cmx1=0.0
|
|
cmy1=0.0
|
|
cmz1=0.0
|
|
cmx2=0.0
|
|
cmy2=0.0
|
|
cmz2=0.0
|
|
summas1=0.0
|
|
summas2=0.0
|
|
do i2=j1,j2
|
|
cmx1=cmx1+c(i2,1)*xmasat(i2)
|
|
cmy1=cmy1+c(i2,2)*xmasat(i2)
|
|
cmz1=cmz1+c(i2,3)*xmasat(i2)
|
|
summas1=summas1+xmasat(i2)
|
|
end do
|
|
cmx1=cmx1/summas1
|
|
cmy1=cmy1/summas1
|
|
cmz1=cmz1/summas1
|
|
if (mdstep.lt.2) then
|
|
rmstrax(i1)=cmx1
|
|
rmstray(i1)=cmy1
|
|
rmstraz(i1)=cmz1
|
|
end if
|
|
if (kdir.le.3) then
|
|
do i2=j3,j4
|
|
cmx2=cmx2+c(i2,1)*xmasat(i2)
|
|
cmy2=cmy2+c(i2,2)*xmasat(i2)
|
|
cmz2=cmz2+c(i2,3)*xmasat(i2)
|
|
summas2=summas2+xmasat(i2)
|
|
end do
|
|
cmx2=cmx2/summas2
|
|
cmy2=cmy2/summas2
|
|
cmz2=cmz2/summas2
|
|
end if
|
|
if (kdir.eq.1) dist=cmx1-cmx2
|
|
if (kdir.eq.2) dist=cmy1-cmy2
|
|
if (kdir.eq.3) dist=cmz1-cmz2
|
|
if (kdir.eq.4) then
|
|
distx=cmx1-rmstrax(i1)
|
|
disty=cmy1-rmstray(i1)
|
|
distz=cmz1-rmstraz(i1)
|
|
dist=sqrt(distx*distx+disty*disty+distz*distz)
|
|
end if
|
|
dismacen(i1)=dist
|
|
dist=dist-rmstra1(i1)
|
|
erh=rmstra2(i1)*dist*dist
|
|
deresdr=2.0*dist*rmstra2(i1)
|
|
* exphu=exp(-rmstra3(i1)*(dist*dist))
|
|
* erh=rmstra2(i1)*(1.0-exphu)
|
|
* deresdr=2.0*rmstra3(i1)*dist*rmstra2(i1)*exphu
|
|
eres=eres+erh
|
|
if (kdir.le.3) then
|
|
do i2=j1,j2
|
|
d(kdir,i2)=d(kdir,i2)+deresdr*xmasat(i2)/summas1
|
|
end do
|
|
do i2=j3,j4
|
|
d(kdir,i2)=d(kdir,i2)-deresdr*xmasat(i2)/summas2
|
|
end do
|
|
end if
|
|
if (kdir.eq.4.and.mdstep.gt.5) then
|
|
do i2=j1,j2
|
|
d(1,i2)=d(1,i2)+deresdr*(distx/dist)*(xmasat(i2)/summas1)
|
|
d(2,i2)=d(2,i2)+deresdr*(disty/dist)*(xmasat(i2)/summas1)
|
|
d(3,i2)=d(3,i2)+deresdr*(distz/dist)*(xmasat(i2)/summas1)
|
|
end do
|
|
end if
|
|
end do
|
|
**********************************************************************
|
|
* *
|
|
* Calculate morphing energy *
|
|
* *
|
|
**********************************************************************
|
|
if (imorph.eq.1) then
|
|
distot=zero
|
|
do i1=1,na
|
|
dmx=c(i1,1)-cmo(i1,1)
|
|
dmy=c(i1,2)-cmo(i1,2)
|
|
dmz=c(i1,3)-cmo(i1,3)
|
|
dism=sqrt(dmx*dmx+dmy*dmy+dmz*dmz)
|
|
distot=distot+dism
|
|
* exphu=exp(-vmo2(i1)*(dism*dism))
|
|
* erh=vmo1(i1)*(1.0-exphu)
|
|
erh=vmo1(i1)*dism
|
|
eres=eres+erh
|
|
* deresddis=2.0*vmo2(i1)*dism*vmo1(i1)*exphu
|
|
deresddis=vmo1(i1)
|
|
drda1=dmx/dism
|
|
drda2=dmy/dism
|
|
drda3=dmz/dism
|
|
d(1,i1)=d(1,i1)+deresddis*drda1
|
|
d(2,i1)=d(2,i1)+deresddis*drda2
|
|
d(3,i1)=d(3,i1)+deresddis*drda3
|
|
end do
|
|
|
|
end if
|
|
|
|
|
|
return
|
|
end
|
|
**********************************************************************
|
|
********************************************************************
|
|
|
|
subroutine calvalres (ja1,ja2,ja3,arg,hr,dhrdc,dargdc)
|
|
|
|
**********************************************************************
|
|
#include "cbka.blk"
|
|
#include "cbkc.blk"
|
|
dimension a(3),b(3),j(3),dradc(3,3),drbdc(3,3),dtdc(3,3),
|
|
$dargdc(3,3),dndc(3,3),dadc(3),dbdc(3),dhrdc(3,3)
|
|
**********************************************************************
|
|
* *
|
|
* Calculate valency angles and their derivatives to cartesian *
|
|
* coordinates for restraint calculations *
|
|
* *
|
|
**********************************************************************
|
|
c$$$* if (ndebug.eq.1) then
|
|
c$$$C* open (65,file='fort.65',status='unknown',access='append')
|
|
c$$$* write (65,*) 'In calvalres'
|
|
c$$$* call timer(65)
|
|
c$$$* close (65)
|
|
c$$$* end if
|
|
|
|
dadc(1)=-1.0
|
|
dadc(2)=1.0
|
|
dadc(3)=0.0
|
|
dbdc(1)=0.0
|
|
dbdc(2)=1.0
|
|
dbdc(3)=-1.0
|
|
do k1=1,3
|
|
do k2=1,3
|
|
dradc(k1,k2)=0.0
|
|
drbdc(k1,k2)=0.0
|
|
end do
|
|
end do
|
|
**********************************************************************
|
|
* *
|
|
* Determine valency angle *
|
|
* *
|
|
**********************************************************************
|
|
call dista2(ja1,ja2,rla,dx1,dy1,dz1)
|
|
call dista2(ja2,ja3,rlb,dx2,dy2,dz2)
|
|
|
|
a(1)=-dx1
|
|
a(2)=-dy1
|
|
a(3)=-dz1
|
|
b(1)=dx2
|
|
b(2)=dy2
|
|
b(3)=dz2
|
|
poem=rla*rlb
|
|
tel=a(1)*b(1)+a(2)*b(2)+a(3)*b(3)
|
|
arg=tel/poem
|
|
arg2=arg*arg
|
|
s1ma22=1.0-arg2
|
|
if (s1ma22.lt.1.0d-10) s1ma22=1.0d-10
|
|
s1ma2=sqrt(s1ma22)
|
|
if (arg.gt.1.0) arg=1.0
|
|
if (arg.lt.-1.0) arg=-1.0
|
|
hr=acos(arg)
|
|
**********************************************************************
|
|
* *
|
|
* Calculate derivative valency angle to cartesian coordinates *
|
|
* *
|
|
**********************************************************************
|
|
do k1=1,3
|
|
dradc(k1,1)=-a(k1)/rla
|
|
dradc(k1,2)=a(k1)/rla
|
|
end do
|
|
|
|
do k1=1,3
|
|
drbdc(k1,2)=b(k1)/rlb
|
|
drbdc(k1,3)=-b(k1)/rlb
|
|
end do
|
|
|
|
do k1=1,3
|
|
do k2=1,3
|
|
dndc(k1,k2)=rla*drbdc(k1,k2)+rlb*dradc(k1,k2)
|
|
dtdc(k1,k2)=a(k1)*dbdc(k2)+b(k1)*dadc(k2)
|
|
dargdc(k1,k2)=(dtdc(k1,k2)-arg*dndc(k1,k2))/poem
|
|
dhrdc(k1,k2)=-dargdc(k1,k2)/s1ma2
|
|
end do
|
|
end do
|
|
|
|
10 continue
|
|
|
|
return
|
|
end
|
|
**********************************************************************
|
|
********************************************************************
|
|
|
|
subroutine calvalhb (ja1,ja2,ja3,ix,iy,iz,arg,hr,dhrdc,dargdc)
|
|
|
|
**********************************************************************
|
|
#include "cbka.blk"
|
|
#include "cbkc.blk"
|
|
dimension a(3),b(3),j(3),dradc(3,3),drbdc(3,3),dtdc(3,3),
|
|
$dargdc(3,3),dndc(3,3),dadc(3),dbdc(3),dhrdc(3,3)
|
|
**********************************************************************
|
|
* *
|
|
* Calculate valency angles and their derivatives to cartesian *
|
|
* coordinates for hydrogen bond calculations *
|
|
* *
|
|
**********************************************************************
|
|
c$$$* if (ndebug.eq.1) then
|
|
c$$$* open (65,file='fort.65',status='unknown',access='append')
|
|
c$$$* write (65,*) 'In calvalhb'
|
|
c$$$* call timer(65)
|
|
c$$$* close (65)
|
|
c$$$* end if
|
|
|
|
dadc(1)=-1.0
|
|
dadc(2)=1.0
|
|
dadc(3)=0.0
|
|
dbdc(1)=0.0
|
|
dbdc(2)=1.0
|
|
dbdc(3)=-1.0
|
|
do k1=1,3
|
|
do k2=1,3
|
|
dradc(k1,k2)=0.0
|
|
drbdc(k1,k2)=0.0
|
|
end do
|
|
end do
|
|
**********************************************************************
|
|
* *
|
|
* Determine valency angle *
|
|
* *
|
|
**********************************************************************
|
|
call dista2(ja1,ja2,rla,dx1,dy1,dz1)
|
|
call dista2(ja2,ja3,rlb,dx2,dy2,dz2)
|
|
|
|
a(1)=-dx1
|
|
a(2)=-dy1
|
|
a(3)=-dz1
|
|
b(1)=dx2
|
|
b(2)=dy2
|
|
b(3)=dz2
|
|
poem=rla*rlb
|
|
tel=a(1)*b(1)+a(2)*b(2)+a(3)*b(3)
|
|
arg=tel/poem
|
|
arg2=arg*arg
|
|
s1ma22=1.0-arg2
|
|
if (s1ma22.lt.1.0d-10) s1ma22=1.0d-10
|
|
s1ma2=sqrt(s1ma22)
|
|
if (arg.gt.1.0) arg=1.0
|
|
if (arg.lt.-1.0) arg=-1.0
|
|
hr=acos(arg)
|
|
**********************************************************************
|
|
* *
|
|
* Calculate derivative valency angle to cartesian coordinates *
|
|
* *
|
|
**********************************************************************
|
|
do k1=1,3
|
|
dradc(k1,1)=-a(k1)/rla
|
|
dradc(k1,2)=a(k1)/rla
|
|
end do
|
|
|
|
do k1=1,3
|
|
drbdc(k1,2)=b(k1)/rlb
|
|
drbdc(k1,3)=-b(k1)/rlb
|
|
end do
|
|
|
|
do k1=1,3
|
|
do k2=1,3
|
|
dndc(k1,k2)=rla*drbdc(k1,k2)+rlb*dradc(k1,k2)
|
|
dtdc(k1,k2)=a(k1)*dbdc(k2)+b(k1)*dadc(k2)
|
|
dargdc(k1,k2)=(dtdc(k1,k2)-arg*dndc(k1,k2))/poem
|
|
dhrdc(k1,k2)=-dargdc(k1,k2)/s1ma2
|
|
end do
|
|
end do
|
|
|
|
10 continue
|
|
|
|
return
|
|
end
|
|
**********************************************************************
|
|
**********************************************************************
|
|
|
|
subroutine caltor(ja1,ja2,ja3,ja4,ht)
|
|
|
|
**********************************************************************
|
|
#include "cbka.blk"
|
|
#include "cbkenergies.blk"
|
|
#include "cbktregime.blk"
|
|
#include "control.blk"
|
|
#include "cbkinit.blk"
|
|
DIMENSION A(3),DRDA(3),DADC(4),DRADC(3,4),DRBDC(3,4),
|
|
$DRCDC(3,4),DHDDC(3,4),DHEDC(3,4),DRVDC(3,4),DTDC(3,4),
|
|
$DNDC(3,4)
|
|
dimension j(4),dvdc1(3,3),dargdc1(3,3),dvdc2(3,3),dargdc2(3,3)
|
|
**********************************************************************
|
|
* *
|
|
* Calculate torsion angle (for internal coordinates output) *
|
|
* *
|
|
**********************************************************************
|
|
c$$$ if (ndebug.eq.1) then
|
|
c$$$C open (65,file='fort.65',status='unknown',access='append')
|
|
c$$$ write (65,*) 'In caltor'
|
|
c$$$ call timer(65)
|
|
c$$$ close (65)
|
|
c$$$ end if
|
|
do k1=1,3
|
|
do k2=1,4
|
|
dhddc(k1,k2)=0.0
|
|
dhedc(k1,k2)=0.0
|
|
dradc(k1,k2)=0.0
|
|
drbdc(k1,k2)=0.0
|
|
drcdc(k1,k2)=0.0
|
|
end do
|
|
end do
|
|
et=0.0
|
|
eco=0.0
|
|
dadc(1)=1.0
|
|
dadc(2)=0.0
|
|
dadc(3)=0.0
|
|
dadc(4)=-1.0
|
|
call dista2(ja1,ja2,rla,dx1,dy1,dz1)
|
|
call dista2(ja2,ja3,rlb,dx2,dy2,dz2)
|
|
call dista2(ja3,ja4,rlc,dx2,dy2,dz2)
|
|
call dista2(ja1,ja4,r4,dx2,dy2,dz2)
|
|
call calvalres(ja1,ja2,ja3,arg1,h1,dvdc1,dargdc1)
|
|
call calvalres(ja2,ja3,ja4,arg2,h2,dvdc2,dargdc2)
|
|
**********************************************************************
|
|
* *
|
|
* Determine torsion angle *
|
|
* *
|
|
**********************************************************************
|
|
d142=r4*r4
|
|
coshd=cos(h1)
|
|
coshe=cos(h2)
|
|
sinhd=sin(h1)
|
|
sinhe=sin(h2)
|
|
poem=2.0*rla*rlc*sinhd*sinhe
|
|
poem2=poem*poem
|
|
tel=rla*rla+rlb*rlb+rlc*rlc-d142-2.0*(rla*rlb*coshd-rla*rlc*
|
|
$coshd*coshe+rlb*rlc*coshe)
|
|
arg=tel/poem
|
|
if (arg.gt.1.0) arg=1.0
|
|
if (arg.lt.-1.0) arg=-1.0
|
|
arg2=arg*arg
|
|
ht=acos(arg)*rdndgr
|
|
k1=ja1
|
|
k2=ja2
|
|
k3=ja3
|
|
k4=ja4
|
|
call dista2(k3,k2,dis,x3,y3,z3)
|
|
y32z32=y3*y3+z3*z3
|
|
wort1=sqrt(y32z32)+1e-6
|
|
wort2=sqrt(y32z32+x3*x3)+1e-6
|
|
sinalf=y3/wort1
|
|
cosalf=z3/wort1
|
|
sinbet=x3/wort2
|
|
cosbet=wort1/wort2
|
|
call dista2(k1,k2,dis,x1,y1,z1)
|
|
x1=x1*cosbet-y1*sinalf*sinbet-z1*cosalf*sinbet
|
|
y1=y1*cosalf-z1*sinalf
|
|
wort3=sqrt(x1*x1+y1*y1)+1e-6
|
|
singam=y1/wort3
|
|
cosgam=x1/wort3
|
|
call dista2(k4,k2,dis,x4,y4,z4)
|
|
x4=x4*cosbet-y4*sinalf*sinbet-z4*cosalf*sinbet
|
|
y4=y4*cosalf-z4*sinalf
|
|
y4=x4*singam-y4*cosgam
|
|
if (y4.gt.0.0) ht=-ht
|
|
if (ht.lt.-179.999999d0) ht=-179.999999d0
|
|
if (ht.gt.179.999999d0) ht=179.999999d0
|
|
|
|
return
|
|
end
|
|
**********************************************************************
|