mirror of
https://github.com/nwchemgit/nwchem.git
synced 2026-07-29 06:35:39 -04:00
New Q-HOP implementation
This commit is contained in:
Tjerk Straatsma
parent
0c22c04e06
commit
2a97e2fc00
16 changed files with 596 additions and 339 deletions
|
|
@ -1,5 +1,5 @@
|
|||
#
|
||||
# $Id: GNUmakefile,v 1.8 2003-07-17 19:44:23 d3j191 Exp $
|
||||
# $Id: GNUmakefile,v 1.9 2006-09-18 15:42:14 d3j191 Exp $
|
||||
#
|
||||
OBJ_OPTIMIZE=cf_d.o cf_fp.o cf_f.o\
|
||||
cf_finish.o\
|
||||
|
|
@ -15,6 +15,7 @@ OBJ_OPTIMIZE=cf_d.o cf_fp.o cf_f.o\
|
|||
cf_shake.o\
|
||||
cf_grid.o\
|
||||
cf_rdtop.o\
|
||||
cf_uhop.o\
|
||||
pme_init.o\
|
||||
pme_fft.o\
|
||||
pme_grid.o
|
||||
|
|
|
|||
|
|
@ -1,4 +1,4 @@
|
|||
c $Id: cf_common.fh,v 1.52 2005-09-14 22:15:56 d3j191 Exp $
|
||||
c $Id: cf_common.fh,v 1.53 2006-09-18 15:42:14 d3j191 Exp $
|
||||
c
|
||||
real*8 zero,tiny,small
|
||||
real*8 fourth,third,half,one,two,three,four,six,twelve
|
||||
|
|
@ -18,6 +18,9 @@ c
|
|||
real*8 avoga,rgas
|
||||
parameter(avoga=6.022045d-07)
|
||||
parameter(rgas=8.31441d-03)
|
||||
c
|
||||
integer mis2
|
||||
parameter(mis2=15)
|
||||
c
|
||||
integer mbs,mpe
|
||||
parameter(mbs=2)
|
||||
|
|
@ -72,7 +75,7 @@ c
|
|||
+ i_s3r1,l_s3r1,i_s3r2,l_s3r2,i_smr1,l_smr1,
|
||||
+ i_smr2,l_smr2,i_smr3,l_smr3,i_smr4,l_smr4,
|
||||
+ i_smr5,l_smr5,i_smr6,l_smr6,
|
||||
+ i_list,l_list,i_snam,l_snam,i_sti,l_sti,
|
||||
+ i_list,l_list,i_snam,l_snam,i_sti,l_sti,i_mprot,l_mprot,
|
||||
+ i_pmew,l_pmew,i_pmes,l_pmes,
|
||||
+ i_theta,l_theta,i_dtheta,l_dtheta,
|
||||
+ i_ipmf,l_ipmf,i_jpmf,l_jpmf,i_rpmf,l_rpmf,i_xpmf,l_xpmf,
|
||||
|
|
@ -96,7 +99,7 @@ c
|
|||
+ i_s3r1,l_s3r1,i_s3r2,l_s3r2,i_smr1,l_smr1,
|
||||
+ i_smr2,l_smr2,i_smr3,l_smr3,i_smr4,l_smr4,
|
||||
+ i_smr5,l_smr5,i_smr6,l_smr6,
|
||||
+ i_list,l_list,i_snam,l_snam,i_sti,l_sti,
|
||||
+ i_list,l_list,i_snam,l_snam,i_sti,l_sti,i_mprot,l_mprot,
|
||||
+ i_pmew,l_pmew,i_pmes,l_pmes,
|
||||
+ i_theta,l_theta,i_dtheta,l_dtheta,
|
||||
+ i_ipmf,l_ipmf,i_jpmf,l_jpmf,i_rpmf,l_rpmf,i_xpmf,l_xpmf,
|
||||
|
|
|
|||
284
src/cafe/cf_f.fh
284
src/cafe/cf_f.fh
|
|
@ -13,7 +13,7 @@ c,rdf,ngt,iagc,jagc,igrc,u,uwmw)
|
|||
c error
|
||||
#endif
|
||||
c
|
||||
c $Id: cf_f.fh,v 1.63 2006-02-07 20:17:04 d3j191 Exp $
|
||||
c $Id: cf_f.fh,v 1.64 2006-09-18 15:42:14 d3j191 Exp $
|
||||
c
|
||||
implicit none
|
||||
c
|
||||
|
|
@ -759,7 +759,7 @@ c
|
|||
c error
|
||||
#endif
|
||||
c
|
||||
c $Id: cf_f.fh,v 1.63 2006-02-07 20:17:04 d3j191 Exp $
|
||||
c $Id: cf_f.fh,v 1.64 2006-09-18 15:42:14 d3j191 Exp $
|
||||
c
|
||||
implicit none
|
||||
c
|
||||
|
|
@ -786,7 +786,7 @@ c
|
|||
integer isal(mscr),isrx(mscr)
|
||||
c
|
||||
integer lswj(*)
|
||||
integer nums
|
||||
integer nums,i
|
||||
integer lswndx(0:msa,2),lswjpt(nums,2),lswin(nums,2)
|
||||
integer list(0:msa)
|
||||
c
|
||||
|
|
@ -821,6 +821,7 @@ c
|
|||
integer is,inum,jnum,jwm,iptr,jwa,iss,isg
|
||||
logical lid,ljd
|
||||
#endif
|
||||
real*8 rtmp
|
||||
c
|
||||
#include "cf_funcs_sfn.fh"
|
||||
#include "bitops_funcs.fh"
|
||||
|
|
@ -833,6 +834,7 @@ c
|
|||
boxi(1)=one/box(1)
|
||||
boxi(2)=one/box(2)
|
||||
boxi(3)=one/box(3)
|
||||
print*,'icntrl=',icntrl
|
||||
if(npbtyp.eq.1.and.nbxtyp.eq.0.and.icntrl.eq.2) then
|
||||
do 101 ipsw=1,npsw
|
||||
do 102 is=1,nums
|
||||
|
|
@ -843,6 +845,7 @@ c
|
|||
inum=lswjpt(is,ipsw)-1
|
||||
if(iand(isss(isa),6).eq.2) iss=-1
|
||||
if(iand(isss(isa),6).eq.4) iss=-1
|
||||
c write(*,'(a,4i5)') 'iss ',isa,isss(isa),iand(isss(isa),6),iss
|
||||
isg=isga(isa)
|
||||
lid=iand(isdt(isa),mdynam).eq.ldynam
|
||||
qi=chg(isq1(isa),1,iset)
|
||||
|
|
@ -1174,6 +1177,9 @@ c isrx(nax+ismn)=isrst
|
|||
c 10 continue
|
||||
c endif
|
||||
c
|
||||
c write(*,'(4i5,2f12.6)')
|
||||
c + lssscl,isga(isa),isss(isfr+isa),iand(isss(isfr+isa),6),
|
||||
c + shift0(1),shift1(1)
|
||||
if(lssscl) then
|
||||
do 10 ismn=1,lswin(isa,ipsw)
|
||||
c isrx(nax+ismn)=isss(isfr+isa)
|
||||
|
|
@ -1285,8 +1291,11 @@ c endif
|
|||
rwx(iax,2)=xi(iax,2)-xj(iax,2,iwa)
|
||||
rwx(iax,3)=xi(iax,3)-xj(iax,3,iwa)
|
||||
rwi2(iax)=one/(rwx(iax,1)**2+rwx(iax,2)**2+rwx(iax,3)**2)
|
||||
rtmp=rwi2(iax)
|
||||
if(isrx(iax).gt.0) rwi2(iax)=one/(one/rwi2(iax)+shift0(1))
|
||||
if(isrx(iax).lt.0) rwi2(iax)=one/(one/rwi2(iax)+shift1(1))
|
||||
c write(*,'(3i5,2f12.6)')
|
||||
c + isga(isal(iax)),isal(iax),isrx(iax),rtmp,rwi2(iax)
|
||||
27 continue
|
||||
c
|
||||
c Lennard-Jones interactions
|
||||
|
|
@ -1889,13 +1898,14 @@ c
|
|||
c error
|
||||
#endif
|
||||
c
|
||||
c $Id: cf_f.fh,v 1.63 2006-02-07 20:17:04 d3j191 Exp $
|
||||
c $Id: cf_f.fh,v 1.64 2006-09-18 15:42:14 d3j191 Exp $
|
||||
c
|
||||
implicit none
|
||||
c
|
||||
#include "cf_common.fh"
|
||||
#include "mafdecls.fh"
|
||||
c
|
||||
real*8 rtmp
|
||||
real*8 xs(msa,3),xsm(msm,3),fs(msa,3,2)
|
||||
real*8 zs(msf,3,3,2),ess(msf,msf,mpe,2)
|
||||
real*8 fss(msf,msf,3,2)
|
||||
|
|
@ -1920,8 +1930,8 @@ c
|
|||
c
|
||||
#if defined(CAFE_FORCES)
|
||||
real*8 dera(6,nsatot)
|
||||
integer lda(10,*)
|
||||
real*8 rda(5,*),uda(4,*)
|
||||
integer lda(16,*)
|
||||
real*8 rda(11,*),uda(4,*)
|
||||
#endif
|
||||
#if defined(CAFE_POLARIZATION)
|
||||
real*8 ps(msa,3,2),psp(msa,3,2,2)
|
||||
|
|
@ -1933,11 +1943,7 @@ c
|
|||
integer ipss,number,isslen,nax,jsaptr
|
||||
integer jnum,lssptr,iax
|
||||
real*8 dercon
|
||||
#if !defined(CAFE_POLARIZATION)
|
||||
integer iss,jss,ilda,ndxqi,ndxqj
|
||||
integer idon,iacc,idong,iaccg,idons,iaccs,idonp,iaccp
|
||||
real*8 uhop(4),xacc(3),racci1
|
||||
#endif
|
||||
c
|
||||
real*8 c6,c12,cf6,cf12
|
||||
real*8 c64,c124
|
||||
real*8 q14,sumen1,sumen2,sumen3
|
||||
|
|
@ -2095,6 +2101,9 @@ c
|
|||
jstt=iand(isss(lssptr),384)
|
||||
if(istt.eq.128.and.jstt.eq.256) isrx(nax+jnum)=999
|
||||
if(istt.eq.256.and.jstt.eq.128) isrx(nax+jnum)=999
|
||||
c
|
||||
c write(*,'(5i5)')
|
||||
c + isga(isfr+isa),isga(lssptr),istt,jstt,isrx(nax+jnum)
|
||||
c
|
||||
endif
|
||||
c
|
||||
|
|
@ -2236,6 +2245,7 @@ c
|
|||
rwi2(iax)=zero
|
||||
isrx(iax)=0
|
||||
endif
|
||||
rtmp=rwi2(iax)
|
||||
if(isrx(iax).gt.0) then
|
||||
rwi2(iax)=one/(one/rwi2(iax)+shift0(1))
|
||||
elseif(isrx(iax).lt.0) then
|
||||
|
|
@ -2244,6 +2254,8 @@ c
|
|||
rwi1(iax)=sqrt(rwi2(iax))
|
||||
isa=isal(iax)
|
||||
jsa=jsal(iax)
|
||||
c write(*,'(3i5,4f12.6)')
|
||||
c + isga(isa),isga(jsa),isrx(iax),shift0(1),shift1(1),rtmp,rwi2(iax)
|
||||
if(jfal(iax).ne.isf) then
|
||||
q=qsa2(isa)*qsa2(jsa)
|
||||
else
|
||||
|
|
@ -2574,177 +2586,6 @@ cx endif
|
|||
nax=nax+lssin(isa,ipss)
|
||||
36 continue
|
||||
c
|
||||
#if defined(CAFE_FORCES)
|
||||
c
|
||||
c qhopping
|
||||
c
|
||||
c lda(1,i) : donor global atom number
|
||||
c lda(2,i) : donor local atom number
|
||||
c lda(3,i) : donor segment number
|
||||
c lda(4,i) : donor protonation state after hopping
|
||||
c lda(5,i) : donor charge mutation
|
||||
c lda(6,i) : acceptor global atom number
|
||||
c lda(7,i) : acceptor local atom number
|
||||
c lda(8,i) : acceptor segment number
|
||||
c lda(9,i) : acceptor protonation state after hopping
|
||||
c lda(10,i) : acceptor charge mutation
|
||||
c
|
||||
if(lhop) then
|
||||
do 2001 ilda=1,nldat
|
||||
c
|
||||
uhop(1)=0.0d0
|
||||
uhop(2)=0.0d0
|
||||
uhop(3)=0.0d0
|
||||
uhop(4)=0.0d0
|
||||
c
|
||||
idon=lda(2,ilda)
|
||||
idong=lda(1,ilda)
|
||||
idons=lda(3,ilda)
|
||||
idonp=lda(4,ilda)
|
||||
iacc=lda(7,ilda)
|
||||
iaccg=lda(6,ilda)
|
||||
iaccs=lda(8,ilda)
|
||||
iaccp=lda(9,ilda)
|
||||
c
|
||||
c write(*,'(a,8i5)') 'DA ',idon,idong,idons,idonp,
|
||||
c + iacc,iaccg,iaccs,iaccp
|
||||
c
|
||||
c Energy of the donor and acceptor segment before hopping
|
||||
c
|
||||
do 2002 iax=1,nax
|
||||
isa=isal(iax)
|
||||
jsa=jsal(iax)
|
||||
iss=isgm(isa)
|
||||
jss=isgm(jsa)
|
||||
c write(*,'(4i4)') isa,jsa,iss,jss
|
||||
c
|
||||
c skip all intrasegment interactions
|
||||
c
|
||||
if(iss.ne.jss.and.(iss.eq.idons.or.jss.eq.idons.or.
|
||||
+ iss.eq.iaccs.or.jss.eq.iaccs)) then
|
||||
c
|
||||
rwi2(iax)=one/(rwx(iax,1)**2+rwx(iax,2)**2+rwx(iax,3)**2)
|
||||
rwi1(iax)=sqrt(rwi2(iax))
|
||||
c write(*,'(4i4,f12.6)') isa,jsa,iss,jss,rwi1(iax)
|
||||
c
|
||||
c donor energy
|
||||
c
|
||||
if((iss.eq.idons.and.jss.ne.iaccs).or.
|
||||
+ (jss.eq.idons.and.iss.ne.iaccs)) then
|
||||
if(jfal(iax).ne.isf) then
|
||||
q=chg(isq2(isa),1,lseq(iss))*
|
||||
+ chg(isq2(jsa),1,lseq(jss))
|
||||
else
|
||||
q=chg(isq3(isa),1,lseq(iss))*
|
||||
+ chg(isq3(jsa),1,lseq(jss))
|
||||
endif
|
||||
uhop(1)=uhop(1)+facu(iax)*q*rwi1(iax)
|
||||
c write(*,'(8i5,2f12.3)')
|
||||
c + isa,lseq(isgm(isa)),isq2(isa),isq3(isa),
|
||||
c + jsa,lseq(isgm(jsa)),isq2(jsa),isq3(jsa),q,uhop(1)
|
||||
endif
|
||||
c
|
||||
c acceptor energy
|
||||
c
|
||||
if((iss.eq.iaccs.and.jss.ne.idons).or.
|
||||
+ (jss.eq.iaccs.and.iss.ne.idons)) then
|
||||
if(jfal(iax).ne.isf) then
|
||||
q=chg(isq2(isa),1,lseq(isgm(isa)))*
|
||||
+ chg(isq2(jsa),1,lseq(isgm(jsa)))
|
||||
else
|
||||
q=chg(isq3(isa),1,lseq(isgm(isa)))*
|
||||
+ chg(isq3(jsa),1,lseq(isgm(jsa)))
|
||||
endif
|
||||
uhop(2)=uhop(2)+facu(iax)*q*rwi1(iax)
|
||||
endif
|
||||
c
|
||||
c donor energy after hopping
|
||||
c
|
||||
if(iss.eq.idons.and.jss.ne.iaccs) then
|
||||
ndxqi=isa
|
||||
if(isga(isa).eq.idong) ndxqi=isa+lda(5,ilda)
|
||||
if(isga(isa).eq.idong+lda(5,ilda)) ndxqi=isa-lda(5,ilda)
|
||||
if(jfal(iax).ne.isf) then
|
||||
q=chg(isq2(ndxqi),1,idonp)*chg(isq2(jsa),1,lseq(isgm(jsa)))
|
||||
else
|
||||
q=chg(isq3(ndxqi),1,idonp)*chg(isq3(jsa),1,lseq(isgm(jsa)))
|
||||
endif
|
||||
uhop(3)=uhop(3)+facu(iax)*q*rwi1(iax)
|
||||
endif
|
||||
if(jss.eq.idons.and.iss.ne.iaccs) then
|
||||
ndxqj=jsa
|
||||
if(isga(jsa).eq.idong) ndxqj=jsa+lda(5,ilda)
|
||||
if(isga(jsa).eq.idong+lda(5,ilda)) ndxqj=jsa-lda(5,ilda)
|
||||
if(jfal(iax).ne.isf) then
|
||||
q=chg(isq2(isa),1,lseq(isgm(isa)))*chg(isq2(ndxqj),1,idonp)
|
||||
else
|
||||
q=chg(isq3(isa),1,lseq(isgm(isa)))*chg(isq3(ndxqj),1,idonp)
|
||||
endif
|
||||
uhop(3)=uhop(3)+facu(iax)*q*rwi1(iax)
|
||||
endif
|
||||
c
|
||||
c acceptor energy after hopping
|
||||
c
|
||||
if(iss.eq.iaccs.and.jss.ne.idons) then
|
||||
if(isga(isa).eq.iaccg) then
|
||||
xacc(1)=rda(1,ilda)-xj(iax,1)
|
||||
xacc(2)=rda(2,ilda)-xj(iax,2)
|
||||
xacc(3)=rda(3,ilda)-xj(iax,3)
|
||||
call cf_pbc(1,xacc,1,xacc,1,0,1,1)
|
||||
else
|
||||
xacc(1)=xi(iax,1)-xj(iax,1)
|
||||
xacc(2)=xi(iax,2)-xj(iax,2)
|
||||
xacc(3)=xi(iax,3)-xj(iax,3)
|
||||
endif
|
||||
racci1=one/sqrt(xacc(1)**2+xacc(2)**2+xacc(3)**2)
|
||||
ndxqi=isa
|
||||
if(isga(isa).eq.iaccg) ndxqi=isa+lda(10,ilda)
|
||||
if(isga(isa).eq.iaccg+lda(10,ilda)) ndxqi=isa-lda(10,ilda)
|
||||
if(jfal(iax).ne.isf) then
|
||||
q=chg(isq2(ndxqi),1,iaccp)*chg(isq2(jsa),1,lseq(isgm(jsa)))
|
||||
else
|
||||
q=chg(isq3(ndxqi),1,iaccp)*chg(isq3(jsa),1,lseq(isgm(jsa)))
|
||||
endif
|
||||
uhop(4)=uhop(4)+facu(iax)*q*racci1
|
||||
endif
|
||||
if(jss.eq.iaccs.and.iss.ne.idons) then
|
||||
if(isga(jsa).eq.iaccg) then
|
||||
xacc(1)=xi(iax,1)-rda(1,ilda)
|
||||
xacc(2)=xi(iax,2)-rda(2,ilda)
|
||||
xacc(3)=xi(iax,3)-rda(3,ilda)
|
||||
call cf_pbc(1,xacc,1,xacc,1,0,1,1)
|
||||
else
|
||||
xacc(1)=xi(iax,1)-xj(iax,1)
|
||||
xacc(2)=xi(iax,2)-xj(iax,2)
|
||||
xacc(3)=xi(iax,3)-xj(iax,3)
|
||||
endif
|
||||
racci1=one/sqrt(xacc(1)**2+xacc(2)**2+xacc(3)**2)
|
||||
ndxqj=jsa
|
||||
if(isga(jsa).eq.iaccg) ndxqj=jsa+lda(5,ilda)
|
||||
if(isga(jsa).eq.iaccg+lda(10,ilda)) ndxqj=jsa-lda(10,ilda)
|
||||
if(jfal(iax).ne.isf) then
|
||||
q=chg(isq2(isa),1,lseq(isgm(isa)))*chg(isq2(ndxqj),1,iaccp)
|
||||
else
|
||||
q=chg(isq3(isa),1,lseq(isgm(isa)))*chg(isq3(ndxqj),1,iaccp)
|
||||
endif
|
||||
uhop(4)=uhop(4)+facu(iax)*q*racci1
|
||||
endif
|
||||
c
|
||||
endif
|
||||
2002 continue
|
||||
c
|
||||
uda(1,ilda)=uda(1,ilda)+uhop(1)
|
||||
uda(2,ilda)=uda(2,ilda)+uhop(2)
|
||||
uda(3,ilda)=uda(3,ilda)+uhop(3)
|
||||
uda(4,ilda)=uda(4,ilda)+uhop(4)
|
||||
c
|
||||
c write(*,'(a,8f12.3)') 'uhop ',(uda(iax,ilda),iax=1,4),uhop
|
||||
c
|
||||
2001 continue
|
||||
endif
|
||||
c
|
||||
#endif
|
||||
c
|
||||
c thermodynamic integration
|
||||
c
|
||||
if(ithint) then
|
||||
|
|
@ -2783,6 +2624,9 @@ c
|
|||
+ rwi2(nax+jnum)*rwi6(nax+jnum)*(c64-c124*rwi6(nax+jnum))
|
||||
endif
|
||||
deriv(15,ipss)=deriv(15,ipss)+dercon
|
||||
c write(*,'(a,3i5,4f12.6)') 'gv ',
|
||||
c + isga(isfr+isa),isga(jsa),isrx(nax+jnum),shift0(4),shift1(4),
|
||||
c + dercon,deriv(15,ipss)
|
||||
#if defined(CAFE_FORCES)
|
||||
if(npgdec.gt.1) then
|
||||
dera(3,isga(isfr+isa))=dera(3,isga(isfr+isa))+half*dercon
|
||||
|
|
@ -2816,6 +2660,9 @@ c
|
|||
if(ireact.ne.0) then
|
||||
deriv(17,ipss)=deriv(17,ipss)+dercon*rffss/rwi2(nax+jnum)
|
||||
endif
|
||||
c write(*,'(a,3i5,4f12.6)') 'gq ',
|
||||
c + isga(isfr+isa),isga(jsa),isrx(nax+jnum),shift0(4),shift1(4),
|
||||
c + dercon,deriv(17,ipss)
|
||||
#if defined(CAFE_FORCES)
|
||||
if(npgdec.gt.1) then
|
||||
dera(4,isga(isfr+isa))=dera(4,isga(isfr+isa))+
|
||||
|
|
@ -3222,7 +3069,7 @@ c
|
|||
+ natoms,ndim,igan,isgm,imol,idyn,ichg,chg,xs,fs,ess,lpbc,lpbcs,
|
||||
+ lupden,lupdti,dera,lseq)
|
||||
c
|
||||
c $Id: cf_f.fh,v 1.63 2006-02-07 20:17:04 d3j191 Exp $
|
||||
c $Id: cf_f.fh,v 1.64 2006-09-18 15:42:14 d3j191 Exp $
|
||||
c
|
||||
c cf_fsb returns forces and energies for solute bonds
|
||||
c
|
||||
|
|
@ -3486,7 +3333,7 @@ c
|
|||
endif
|
||||
c
|
||||
if(ipme.ne.0) then
|
||||
qij=chg(ichg(isa),1,1)*chg(ichg(jsa),1,1)
|
||||
qij=chg(ichg(isa),1,ibset)*chg(ichg(jsa),1,ibset)
|
||||
ferfc=one-erfc(ealpha*rss)
|
||||
fderfc=-(ealpha*derfc(ealpha*rss))
|
||||
epmecs=epmecs-ferfc*qij*rssi
|
||||
|
|
@ -3532,7 +3379,7 @@ c
|
|||
+ natoms,ndim,igan,isgm,imol,idyn,ichg,chg,xs,fs,ess,lpbc,lpbcs,
|
||||
+ lupden,lupdti,dera,lseq)
|
||||
c
|
||||
c $Id: cf_f.fh,v 1.63 2006-02-07 20:17:04 d3j191 Exp $
|
||||
c $Id: cf_f.fh,v 1.64 2006-09-18 15:42:14 d3j191 Exp $
|
||||
c
|
||||
c cf_fsh returns forces and energies for solute angles
|
||||
c
|
||||
|
|
@ -3955,7 +3802,7 @@ c
|
|||
+ natoms,ndim,igan,isgm,imol,idyn,xs,fs,ess,lpbc,lpbcs,lupden,
|
||||
+ lupdti,dera,lseq)
|
||||
c
|
||||
c $Id: cf_f.fh,v 1.63 2006-02-07 20:17:04 d3j191 Exp $
|
||||
c $Id: cf_f.fh,v 1.64 2006-09-18 15:42:14 d3j191 Exp $
|
||||
c
|
||||
c cf_fsd returns forces and energies for solute angles
|
||||
c
|
||||
|
|
@ -4301,7 +4148,7 @@ c
|
|||
+ natoms,ndim,igan,isgm,imol,idyn,xs,fs,ess,lpbc,lpbcs,lupden,
|
||||
+ lupdti,dera,lseq)
|
||||
c
|
||||
c $Id: cf_f.fh,v 1.63 2006-02-07 20:17:04 d3j191 Exp $
|
||||
c $Id: cf_f.fh,v 1.64 2006-09-18 15:42:14 d3j191 Exp $
|
||||
c
|
||||
c cf_fso returns forces and energies for solute improper dihedral angles
|
||||
c
|
||||
|
|
@ -4628,9 +4475,10 @@ c
|
|||
return
|
||||
end
|
||||
subroutine cf_fst(nthrds,indexl,mst,idtn,vdw,chg,natoms,ndim,
|
||||
+ iatt,igan,isgm,imol,idyn,ichg,xs,fs,ess,lpbc,lpbcs,dera,esa,lseq)
|
||||
+ iatt,igan,isgm,imol,idyn,ichg,isss,xs,fs,ess,lpbc,lpbcs,dera,
|
||||
+ esa,lseq)
|
||||
c
|
||||
c $Id: cf_f.fh,v 1.63 2006-02-07 20:17:04 d3j191 Exp $
|
||||
c $Id: cf_f.fh,v 1.64 2006-09-18 15:42:14 d3j191 Exp $
|
||||
c
|
||||
c cf_fst returns forces and energies for solute third neighbors
|
||||
c
|
||||
|
|
@ -4682,6 +4530,7 @@ c
|
|||
integer natoms,ndim
|
||||
real*8 vdw(mat,mat,map,mset),chg(mqt,mqp,mset)
|
||||
integer iatt(ndim),igan(ndim),imol(ndim),idyn(ndim),ichg(ndim)
|
||||
integer isss(msa)
|
||||
real*8 xs(ndim,3),fs(ndim,3)
|
||||
real*8 ess(msf,msf,mpe,2),esa(nsa)
|
||||
c
|
||||
|
|
@ -4694,7 +4543,7 @@ c
|
|||
c declaration of local variables
|
||||
c ------------------------------
|
||||
c
|
||||
integer i,j,ist,isa,jsa,isf,jsf,isad
|
||||
integer i,j,ist,isa,jsa,isf,jsf,isad,istt,jstt,isrx
|
||||
real*8 factu,etermq,eterml
|
||||
real*8 rxx,rxy,rxz,r2,r2i,r1i,dfs,dfsu,r6i,c6,c12,cf6,cf12,q
|
||||
real*8 ferfc,fderfc,qij,dfsp,c64,c124,dercon,qi,qi4,qj,qj4,q14
|
||||
|
|
@ -4746,6 +4595,30 @@ c find solute molecule numbers involved
|
|||
c
|
||||
isf=imol(isa)
|
||||
jsf=imol(jsa)
|
||||
c
|
||||
if(lssscl) then
|
||||
c
|
||||
istt=iand(isss(isa),48)
|
||||
jstt=iand(isss(jsa),48)
|
||||
if(isf.ne.jsf) then
|
||||
if(istt.eq.16.or.jstt.eq.16) isrx=-1
|
||||
if(istt.eq.32.or.jstt.eq.32) isrx=1
|
||||
endif
|
||||
c
|
||||
istt=iand(isss(isa),384)
|
||||
jstt=iand(isss(jsa),384)
|
||||
if(istt.eq.128.or.jstt.eq.128) isrx=-2
|
||||
if(istt.eq.256.or.jstt.eq.256) isrx=2
|
||||
c
|
||||
istt=iand(isss(isa),384)
|
||||
jstt=iand(isss(jsa),384)
|
||||
if(istt.eq.128.and.jstt.eq.256) isrx=999
|
||||
if(istt.eq.256.and.jstt.eq.128) isrx=999
|
||||
c
|
||||
c write(*,'(5i5)')
|
||||
c + isga(isfr+isa),isga(lssptr),istt,jstt,isrx(nax+jnum)
|
||||
c
|
||||
endif
|
||||
c
|
||||
c determine energy factor depending on atoms being dynamic
|
||||
c or fixed
|
||||
|
|
@ -4779,6 +4652,19 @@ c
|
|||
c
|
||||
r2=rxx*rxx+rxy*rxy+rxz*rxz
|
||||
r2i=one/r2
|
||||
c
|
||||
if(lssscl) then
|
||||
if(isrx.eq.999) then
|
||||
r2i=zero
|
||||
isrx=0
|
||||
endif
|
||||
if(isrx.gt.0) then
|
||||
r2i=one/(one/r2i+shift0(1))
|
||||
elseif(isrx.lt.0) then
|
||||
r2i=one/(one/r2i+shift1(1))
|
||||
endif
|
||||
endif
|
||||
c
|
||||
r1i=sqrt(r2i)
|
||||
dfs=zero
|
||||
dfsu=zero
|
||||
|
|
@ -4894,6 +4780,15 @@ c c124=cb12(iatt(isa),iatt(jsa),4)
|
|||
c64=vdw(iatt(isa),iatt(jsa),2,4)
|
||||
c124=vdw(iatt(isa),iatt(jsa),4,4)
|
||||
dercon=(c124*r6i-c64)*r6i
|
||||
if(isrx.gt.0) then
|
||||
c64=three*vdw(iatt(isa),iatt(jsa),1,iset)
|
||||
c124=six*vdw(iatt(isa),iatt(jsa),3,iset)
|
||||
dercon=dercon+shift0(4)*r2i*r6i*(c64-c124*r6i)
|
||||
elseif(isrx.lt.0) then
|
||||
c64=three*vdw(iatt(isa),iatt(jsa),1,iset)
|
||||
c124=six*vdw(iatt(isa),iatt(jsa),3,iset)
|
||||
dercon=dercon+shift1(4)*r2i*r6i*(c64-c124*r6i)
|
||||
endif
|
||||
deriv(15,1)=deriv(15,1)+dercon
|
||||
#if defined(CAFE_FORCES)
|
||||
if(npgdec.gt.1) then
|
||||
|
|
@ -4915,6 +4810,13 @@ c qj4=qs(jsa,4)
|
|||
dercon=q14fac*(qi*qj4+qj*qi4)*r1i
|
||||
else
|
||||
dercon=q14fac*(qi*qj4+qj*qi4)*r1i
|
||||
endif
|
||||
if(isrx.gt.1) then
|
||||
dercon=dercon-half*shift0(4)*
|
||||
+ chg(ichg(isa),1,iset)*chg(ichg(jsa),1,iset)*r2i
|
||||
elseif(isrx.lt.-1) then
|
||||
dercon=dercon-half*shift1(4)*
|
||||
+ chg(ichg(isa),1,iset)*chg(ichg(jsa),1,iset)*r2i
|
||||
endif
|
||||
deriv(17,1)=deriv(17,1)+dercon
|
||||
if(ireact.ne.0) then
|
||||
|
|
@ -4993,7 +4895,7 @@ c
|
|||
+ mwd,nwd,ndp,idih,dih,rdih,mwo,nwo,nop,iimp,dimp,rimp,
|
||||
+ mwt,nwt,idwt,mwn,nwn,idwn)
|
||||
c
|
||||
c $Id: cf_f.fh,v 1.63 2006-02-07 20:17:04 d3j191 Exp $
|
||||
c $Id: cf_f.fh,v 1.64 2006-09-18 15:42:14 d3j191 Exp $
|
||||
c
|
||||
implicit none
|
||||
c
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
subroutine cf_finish()
|
||||
c $Id: cf_finish.F,v 1.22 2003-11-25 16:40:12 d3j191 Exp $
|
||||
c $Id: cf_finish.F,v 1.23 2006-09-18 15:42:14 d3j191 Exp $
|
||||
implicit none
|
||||
c
|
||||
#include "cf_common.fh"
|
||||
|
|
@ -142,6 +142,11 @@ c
|
|||
if(.not.ma_pop_stack(l_sti))
|
||||
+ call md_abort('Failed to deallocate sti',0)
|
||||
c
|
||||
if(lqhop) then
|
||||
if(.not.ma_pop_stack(l_mprot))
|
||||
+ call md_abort('Failed to deallocate mprot',0)
|
||||
endif
|
||||
c
|
||||
if(.not.ma_pop_stack(l_snam))
|
||||
+ call md_abort('Failed to deallocate snam',0)
|
||||
c
|
||||
|
|
|
|||
|
|
@ -1,4 +1,4 @@
|
|||
c $Id: cf_forces.fh,v 1.23 2005-01-08 19:17:20 d3j191 Exp $
|
||||
c $Id: cf_forces.fh,v 1.24 2006-09-18 15:42:14 d3j191 Exp $
|
||||
#if defined(CAFE_POLARIZATION)
|
||||
subroutine forcep(lself,lpbcs,
|
||||
+ xw,xwm,fw,pw,pwp,zwi,rtos,iwdt,iwz,
|
||||
|
|
@ -349,7 +349,7 @@ c
|
|||
if(npairs.gt.0) then
|
||||
call cf_fst(npairs,int_mb(lptr+1),mtt(2),int_mb(i_itrd(2)),
|
||||
+ dbl_mb(i_vdw),dbl_mb(i_chg),
|
||||
+ max(isto,jsto),msa,isat,isga,isgm,ismf,isdt,isq3,
|
||||
+ max(isto,jsto),msa,isat,isga,isgm,ismf,isdt,isq3,isss,
|
||||
+ xs,fs,ess,lpbc,lpbcs,dbl_mb(i_sti),esa,int_mb(i_lseq))
|
||||
lptr=lptr+1+npairs
|
||||
else
|
||||
|
|
|
|||
|
|
@ -1,7 +1,7 @@
|
|||
logical function cf_hopping(stimei,issgm,isgan,ishop,
|
||||
+ xs,nsaloc)
|
||||
logical function cf_hopping(lpbc,lpbcs,stimei,
|
||||
+ isl,issgm,isgan,isq3,ishop,xs,nsaloc)
|
||||
c
|
||||
c $Id: cf_hopping.F,v 1.16 2005-10-03 21:03:12 edo Exp $
|
||||
c $Id: cf_hopping.F,v 1.17 2006-09-18 15:42:14 d3j191 Exp $
|
||||
c
|
||||
implicit none
|
||||
c
|
||||
|
|
@ -13,13 +13,16 @@ c
|
|||
logical cf_hop
|
||||
external cf_hop
|
||||
c
|
||||
integer issgm(msa),nsaloc,isgan(msa),ishop(msa)
|
||||
logical lpbc,lpbcs
|
||||
integer isl(msa,mis2)
|
||||
integer issgm(msa),nsaloc,isgan(msa),ishop(msa),isq3(msa)
|
||||
real*8 xs(msa,3)
|
||||
real*8 stimei
|
||||
c
|
||||
integer i_itmp,l_itmp,i_dtmp,l_dtmp
|
||||
c
|
||||
stime=stimei
|
||||
cf_hopping=.false.
|
||||
c
|
||||
if(lpair.and.nhop.gt.0) then
|
||||
lhop=.true.
|
||||
|
|
@ -27,13 +30,20 @@ c
|
|||
if(.not.ma_push_get(mt_int,np,'itmp',l_itmp,i_itmp))
|
||||
+ call md_abort('Failed to allocate itmp',me)
|
||||
c
|
||||
call cf_hoplist(int_mb(i_itmp),int_mb(i_lda),dbl_mb(i_rda))
|
||||
call cf_hoplist(issgm,nsaloc,int_mb(i_itmp),int_mb(i_lda),
|
||||
+ dbl_mb(i_rda))
|
||||
c
|
||||
if(.not.ma_pop_stack(l_itmp))
|
||||
+ call md_abort('Failed to deallocate itmp',me)
|
||||
c
|
||||
cf_hopping=.true.
|
||||
elseif(lhop) then
|
||||
c
|
||||
c evaluate the E12
|
||||
c
|
||||
call cf_uhop(lpbc,lpbcs,int_mb(i_lseq),isl,issgm,isgan,isq3,
|
||||
+ xs,nsaloc,
|
||||
+ dbl_mb(i_chg),int_mb(i_lda),dbl_mb(i_rda),dbl_mb(i_uda))
|
||||
c
|
||||
if(.not.ma_push_get(mt_dbl,nldat,'dtmp',l_dtmp,i_dtmp))
|
||||
+ call md_abort('Failed to allocate dtmp',me)
|
||||
|
|
@ -54,7 +64,7 @@ c
|
|||
c
|
||||
return
|
||||
end
|
||||
subroutine cf_hoplist(nhopl,lda,rda)
|
||||
subroutine cf_hoplist(issgm,nsaloc,nhopl,lda,rda)
|
||||
c
|
||||
implicit none
|
||||
c
|
||||
|
|
@ -62,10 +72,11 @@ c
|
|||
#include "global.fh"
|
||||
#include "mafdecls.fh"
|
||||
#include "msgids.fh"
|
||||
#include "util.fh"
|
||||
c
|
||||
integer nhopl(np)
|
||||
integer lda(10,*)
|
||||
real*8 rda(5,*)
|
||||
integer issgm(msa),nsaloc,nhopl(np)
|
||||
integer lda(16,*)
|
||||
real*8 rda(11,*)
|
||||
c
|
||||
integer i,j,ioff
|
||||
c
|
||||
|
|
@ -73,52 +84,82 @@ c
|
|||
nhopl(i)=0
|
||||
1 continue
|
||||
nhopl(me+1)=nlda
|
||||
c write(*,'(i5,a,i5)') me,' nlda= ',nlda
|
||||
c
|
||||
call ga_igop(mcf_76,nhopl,np,'+')
|
||||
c
|
||||
do 2 i=2,np
|
||||
nhopl(i)=nhopl(i)+nhopl(i-1)
|
||||
2 continue
|
||||
c if(me.eq.0) write(*,'(a,10i5)') 'nhopl ',(nhopl(i),i=1,np)
|
||||
c write(*,'(i5,a,10i5)') me,' nhopl ',(nhopl(i),i=1,np)
|
||||
c
|
||||
c do 222 i=1,nlda
|
||||
c write(*,'(12i5)') me,i,(lda(j,i),j=1,10)
|
||||
c 222 continue
|
||||
c
|
||||
ioff=0
|
||||
if(me.gt.0) then
|
||||
ioff=nhopl(me)
|
||||
do 3 i=1,nlda
|
||||
do 4 j=1,10
|
||||
do 4 j=1,16
|
||||
lda(j,i+ioff)=lda(j,i)
|
||||
4 continue
|
||||
do 5 j=1,5
|
||||
do 5 j=1,8
|
||||
rda(j,i+ioff)=rda(j,i)
|
||||
5 continue
|
||||
3 continue
|
||||
do 6 i=1,nhopl(me)
|
||||
do 7 j=1,10
|
||||
do 7 j=1,16
|
||||
lda(j,i)=0
|
||||
7 continue
|
||||
do 8 j=1,5
|
||||
do 8 j=1,8
|
||||
rda(j,i)=0.0d0
|
||||
8 continue
|
||||
6 continue
|
||||
endif
|
||||
if(me.lt.np-1) then
|
||||
do 9 i=nhopl(me+1)+1,nhopl(np)
|
||||
do 10 j=1,10
|
||||
do 10 j=1,16
|
||||
lda(j,i)=0
|
||||
10 continue
|
||||
do 11 j=1,5
|
||||
do 11 j=1,8
|
||||
rda(j,i)=0.0d0
|
||||
11 continue
|
||||
9 continue
|
||||
endif
|
||||
c
|
||||
nldat=nhopl(np)
|
||||
call ga_igop(mcf_77,lda,10*nldat,'+')
|
||||
call ga_dgop(mcf_78,rda,5*nldat,'+')
|
||||
call ga_igop(mcf_77,lda,16*nldat,'+')
|
||||
c
|
||||
c do 322 i=1,nldat
|
||||
c write(*,'(18i5)') me,i,(lda(j,i),j=1,16)
|
||||
c 322 continue
|
||||
c
|
||||
do 201 i=1,nldat
|
||||
do 202 j=1,nsaloc
|
||||
if(lda(3,i).eq.issgm(j)) then
|
||||
if(lda(12,i).eq.0) lda(12,i)=j
|
||||
lda(13,i)=j
|
||||
endif
|
||||
if(lda(8,i).eq.issgm(j)) then
|
||||
if(lda(15,i).eq.0) lda(15,i)=j
|
||||
lda(16,i)=j
|
||||
endif
|
||||
202 continue
|
||||
201 continue
|
||||
c
|
||||
if(me.gt.0) then
|
||||
do 221 j=1,nldat
|
||||
do 222 i=1,11
|
||||
lda(i,j)=0
|
||||
222 continue
|
||||
lda(14,j)=0
|
||||
221 continue
|
||||
endif
|
||||
c
|
||||
call ga_igop(mcf_77,lda,16*nldat,'+')
|
||||
call ga_dgop(mcf_78,rda,11*nldat,'+')
|
||||
c
|
||||
c do 2222 i=1,nlda
|
||||
c write(*,'(18i5)') me,i,(lda(j,i),j=1,16)
|
||||
c 2222 continue
|
||||
c
|
||||
nlda=0
|
||||
do 12 i=1,nldat
|
||||
|
|
@ -126,24 +167,31 @@ c
|
|||
if(lda(1,i).eq.lda(1,j).and.lda(6,i).eq.lda(6,j)) goto 12
|
||||
13 continue
|
||||
nlda=nlda+1
|
||||
do 14 j=1,10
|
||||
do 14 j=1,16
|
||||
lda(j,nlda)=lda(j,i)
|
||||
14 continue
|
||||
do 15 j=1,5
|
||||
do 15 j=1,8
|
||||
rda(j,nlda)=rda(j,i)
|
||||
15 continue
|
||||
12 continue
|
||||
nldat=nlda
|
||||
c
|
||||
ctest
|
||||
c if(me.eq.0) then
|
||||
c print*,'Number of donor-acceptor pairs is ',nldat
|
||||
c do 12 i=1,nldat
|
||||
c write(*,111) (lda(j,i),j=1,10),(rda(j,i),j=1,3)
|
||||
c 111 format(10x,10i5,3f12.6)
|
||||
c 12 continue
|
||||
c endif
|
||||
ctest
|
||||
c
|
||||
if(me.eq.0) then
|
||||
if(util_print('qhop',print_high)) then
|
||||
write(lfnhop,110) nldat
|
||||
110 format(/,'Number of donor-acceptor pairs is ',i5,/,
|
||||
+ '---------- Donor --------- --------- Acceptor ----',
|
||||
+ ' --- Donor --- ---Acceptor---',/,
|
||||
+ ' glob loc sgm prot off glob loc sgm prot off',
|
||||
+ ' own from to own from to',/)
|
||||
do 212 i=1,nldat
|
||||
write(lfnhop,111) (lda(j,i),j=1,16),(rda(j,i),j=1,3)
|
||||
111 format(16i5,3f12.6)
|
||||
212 continue
|
||||
endif
|
||||
endif
|
||||
c
|
||||
return
|
||||
end
|
||||
logical function cf_hop(lseq,issgm,lda,rda,uda,ndx,pda,
|
||||
|
|
@ -160,40 +208,40 @@ c
|
|||
#include "bitops_funcs.fh"
|
||||
c
|
||||
integer lseq(mseq),issgm(msa)
|
||||
integer lda(10,*),nsaloc,isgan(msa),ishop(msa)
|
||||
real*8 rda(5,*),uda(4,*),pda(*),xs(msa,3),x(3)
|
||||
integer lda(16,*),nsaloc,isgan(msa),ishop(msa)
|
||||
real*8 rda(11,*),uda(4,*),pda(*),xs(msa,3),x(3)
|
||||
integer ndx(*)
|
||||
integer lsthop(2,*)
|
||||
real*8 timhop(*)
|
||||
c
|
||||
integer i,j,k,ndxrel
|
||||
real*8 dprob,drand,e12
|
||||
real*8 dprob,drand,e12,e120
|
||||
c
|
||||
call ga_dgop(mcf_79,uda,4*nldat,'+')
|
||||
c
|
||||
if(me.eq.0) then
|
||||
if(util_print('qhop',print_high)) then
|
||||
write(lfnhop,110)
|
||||
110 format(/,
|
||||
write(lfnhop,110) stime
|
||||
110 format(/,' Time ',f12.6,/,
|
||||
+ ' QHOP ----donor--- --acceptor-- ',
|
||||
+ ' pre-hop energies post-hop energies ',/,
|
||||
+ ' numb atom res prot atom res prot',
|
||||
+ ' donor acceptor donor acceptor ',
|
||||
+ ' e12 distance angle probability',/)
|
||||
+ ' e12 e120 distance angle probability',/)
|
||||
endif
|
||||
do 1 i=1,nldat
|
||||
c write(*,'(i5,5f12.6)') i,(rda(j,i),j=1,5)
|
||||
c write(*,'(i5,8f12.6)') i,(rda(j,i),j=1,8)
|
||||
call qhop_prob(lda(1,i),lda(3,i),lseq(lda(3,i)),
|
||||
+ lda(6,i),lda(8,i),lseq(lda(8,i)),uda(1,i),rda(4,i),pda(i))
|
||||
c write(*,'(i5,5f12.6)') i,(rda(j,i),j=1,5)
|
||||
+ lda(6,i),lda(8,i),lseq(lda(8,i)),uda(1,i),rda(4,i),pda(i),e120)
|
||||
c write(*,'(i5,8f12.6)') i,(rda(j,i),j=1,8)
|
||||
c pda(i)=0.0d0
|
||||
c
|
||||
if(util_print('qhop',print_high)) then
|
||||
e12=uda(3,i)+uda(4,i)-uda(1,i)-uda(2,i)
|
||||
write(lfnhop,111) i,lda(1,i),lda(3,i),lseq(lda(3,i)),
|
||||
+ lda(6,i),lda(8,i),lseq(lda(8,i)),
|
||||
+ (uda(j,i),j=1,4),e12,rda(4,i),rda(5,i),pda(i)
|
||||
111 format(7i5,7f12.6,f12.9)
|
||||
+ (uda(j,i),j=1,4),e12,e120*4.184,rda(4,i),rda(5,i),pda(i)
|
||||
111 format(7i5,8f12.6,f12.9)
|
||||
endif
|
||||
c
|
||||
1 continue
|
||||
|
|
|
|||
|
|
@ -7,7 +7,7 @@
|
|||
+ ipset2,issscl,delta,nfanal,lpbc,npgi,fldi,fvec,ffrq,npenrg,ictrl,
|
||||
+ nbiasi,mropti,incl,ltwn,nseqi,i_lseqi,nfhopi,rhopi,thopi,ndumsi,
|
||||
+ ipbtpi,lfnhopi,iradgi,nbgeti,npreci)
|
||||
c $Id: cf_init.F,v 1.108 2005-09-14 22:15:56 d3j191 Exp $
|
||||
c $Id: cf_init.F,v 1.109 2006-09-18 15:42:14 d3j191 Exp $
|
||||
implicit none
|
||||
c
|
||||
#include "cf_common.fh"
|
||||
|
|
@ -228,6 +228,11 @@ c
|
|||
c
|
||||
if(.not.ma_push_get(mt_byte,16*nsatot,'snam',l_snam,i_snam))
|
||||
+ call md_abort('Failed to allocate snam',me)
|
||||
c
|
||||
if(lqhop) then
|
||||
if(.not.ma_push_get(mt_int,mseq,'mprot',l_mprot,i_mprot))
|
||||
+ call md_abort('Failed to allocate mprot',me)
|
||||
endif
|
||||
c
|
||||
if(npgdec.gt.1) then
|
||||
if(.not.ma_push_get(mt_dbl,6*nsatot,'sti',l_sti,i_sti))
|
||||
|
|
@ -278,9 +283,9 @@ c proton hopping: donor-acceptor pair list allocation
|
|||
c ---------------------------------------------------
|
||||
c
|
||||
if(nhop.gt.0) then
|
||||
if(.not.ma_push_get(mt_int,10*nhop*3,'lda',l_lda,i_lda))
|
||||
if(.not.ma_push_get(mt_int,16*nhop*3,'lda',l_lda,i_lda))
|
||||
+ call md_abort('Failed to allocate lda',0)
|
||||
if(.not.ma_push_get(mt_dbl,5*nhop*3,'rda',l_rda,i_rda))
|
||||
if(.not.ma_push_get(mt_dbl,11*nhop*3,'rda',l_rda,i_rda))
|
||||
+ call md_abort('Failed to allocate rda',0)
|
||||
if(.not.ma_push_get(mt_dbl,4*nhop*3,'uda',l_uda,i_uda))
|
||||
+ call md_abort('Failed to allocate uda',0)
|
||||
|
|
@ -704,7 +709,7 @@ c
|
|||
#include "mafdecls.fh"
|
||||
c
|
||||
integer ia
|
||||
real*8 chgi(24,2)
|
||||
real*8 chgi(24,3)
|
||||
c
|
||||
integer i
|
||||
c
|
||||
|
|
@ -730,7 +735,7 @@ c
|
|||
#include "mafdecls.fh"
|
||||
c
|
||||
integer ia
|
||||
real*8 chgi(24,2),chg(mqt,mqp,mset)
|
||||
real*8 chgi(24,3),chg(mqt,mqp,mset)
|
||||
c
|
||||
integer i
|
||||
c
|
||||
|
|
@ -738,11 +743,13 @@ c
|
|||
do 1 i=1,nparms
|
||||
chg(ia,1,i)=qfac*chgi(i,1)
|
||||
chg(ia,2,i)=chgi(i,2)
|
||||
chg(ia,3,i)=qfac*chgi(i,3)
|
||||
1 continue
|
||||
else
|
||||
do 2 i=1,3
|
||||
chg(ia,1,i)=qfac*chgi(i,1)
|
||||
chg(ia,2,i)=chgi(i,2)
|
||||
chg(ia,3,i)=qfac*chgi(i,3)
|
||||
2 continue
|
||||
do 3 i=1,2
|
||||
chg(ia,i,4)=chg(ia,i,3)-chg(ia,i,2)
|
||||
|
|
@ -1082,6 +1089,19 @@ c
|
|||
call cf_pcopy(idp,ip,p,6,mdp(id),int_mb(i_idih(id)),
|
||||
+ dbl_mb(i_dih(id)),mdt(id),
|
||||
+ ip2(9+(id-1)*12),ip3(9+(id-1)*12),ith(9+(id-1)*12))
|
||||
c
|
||||
return
|
||||
end
|
||||
subroutine cf_parseq(id,idp)
|
||||
c
|
||||
c cf_parseq
|
||||
c
|
||||
#include "cf_common.fh"
|
||||
#include "mafdecls.fh"
|
||||
c
|
||||
integer id,idp
|
||||
c
|
||||
int_mb(i_mprot+id-1)=idp
|
||||
c
|
||||
return
|
||||
end
|
||||
|
|
|
|||
|
|
@ -1,6 +1,6 @@
|
|||
subroutine cf_lww(lself,lpbc,xwm,idt,iwfr,iwto,jwfr,jwto,
|
||||
+ nlocw,mpairs,npairs,lwwjpt,lwwin,lwwj,list,rwx,rw)
|
||||
c $Id: cf_list.F,v 1.32 2005-01-08 19:17:20 d3j191 Exp $
|
||||
c $Id: cf_list.F,v 1.33 2006-09-18 15:42:14 d3j191 Exp $
|
||||
implicit none
|
||||
c
|
||||
#include "cf_common.fh"
|
||||
|
|
@ -351,7 +351,7 @@ c
|
|||
subroutine cf_lss(lself,lpbc,lpbcs,xs,iga,idt,igr,isgm,iss,ish,
|
||||
+ isz,
|
||||
+ isfr,isto,jsfr,jsto,nlocs,mpairs,npairs,lssjpt,lssin,lssj,list,
|
||||
+ rwx,rw,ndxt,nt,ndxx,nx,lda,rda,lseq,lsthop)
|
||||
+ rwx,rw,ndxt,nt,ndxx,nx,lda,rda,lseq,lsthop,mprot)
|
||||
c
|
||||
implicit none
|
||||
c
|
||||
|
|
@ -368,19 +368,19 @@ c
|
|||
integer lssj(*),lssjpt(nlocs,2),lssin(nlocs,2)
|
||||
integer list(msa,3),iga(msa),igr(msa),idt(msa),iss(msa),isz(msa)
|
||||
integer ish(msa),isgm(msa)
|
||||
integer lseq(mseq),lsthop(2,*)
|
||||
logical lself,lpbc,lpbcs
|
||||
integer lseq(mseq),lsthop(2,*),mprot(mseq)
|
||||
logical lself,lpbc,lpbcs,lsame
|
||||
c
|
||||
integer isafr,isfr,isato,isto,jsafr,jsfr,jsato,jsto
|
||||
integer mpairs,npairs,isa,jsa,nlist1,nlist2,idoff,iaoff
|
||||
integer ifrsx,itosx,ifrst,itost,jfrsx,jtosx,jfrst,jtost
|
||||
integer iskip1,iskip2,istsum,ist,jst,isxsum,isx,jsx,ilist
|
||||
integer nlist0,isa0,iisa,isglob,ipss,jsglob,nlist,leq,nssl0
|
||||
integer lda(10,*)
|
||||
real*8 rda(5,*)
|
||||
integer lda(16,*)
|
||||
real*8 rda(11,*)
|
||||
c
|
||||
real*8 xacc(5)
|
||||
integer idon,iacc,i
|
||||
real*8 xacc(11)
|
||||
integer idon,iacc,i,idown,idfr,idto,iaown,iafr,iato
|
||||
c
|
||||
c
|
||||
c this subroutine evaluates the solute-solute pairlist
|
||||
|
|
@ -391,6 +391,8 @@ c
|
|||
jsato=jsto
|
||||
if(lself) isato=isato-1
|
||||
c
|
||||
c write(*,'(i5,a,i5)') me,' cf_lss nseq=',nseq
|
||||
c write(*,'(20i5)') (lseq(i),i=1,nseq)
|
||||
c
|
||||
npairs=0
|
||||
c
|
||||
|
|
@ -570,11 +572,15 @@ c
|
|||
if(nhop.gt.0) then
|
||||
idon=0
|
||||
iacc=0
|
||||
idown=0
|
||||
iaown=0
|
||||
if(ish(isa).gt.0) then
|
||||
if(iand(ish(isa),1).eq.0) then
|
||||
idon=isa
|
||||
idown=me
|
||||
else
|
||||
iacc=isa
|
||||
iaown=me
|
||||
endif
|
||||
endif
|
||||
if(idon.gt.0.or.iacc.gt.0) then
|
||||
|
|
@ -585,29 +591,39 @@ c
|
|||
if(ish(jsa).gt.0) then
|
||||
if(iacc.gt.0.and.iand(ish(jsa),1).eq.0) then
|
||||
idon=jsa
|
||||
if(lself) idown=me
|
||||
elseif(idon.gt.0.and.iand(ish(jsa),1).eq.1) then
|
||||
iacc=jsa
|
||||
if(lself) iaown=me
|
||||
endif
|
||||
endif
|
||||
if(idon.gt.0.and.iacc.gt.0.and.rw(jsa).lt.rhop2) then
|
||||
ctmp
|
||||
c temporary change to exclude same to same protonation state
|
||||
c
|
||||
if(lseq(isgm(idon)).ne.lseq(isgm(iacc))) then
|
||||
ctmp
|
||||
if(lseq(isgm(idon))+ish(idon)/2.le.mprot(idon)
|
||||
+ .or.lseq(isgm(iacc))-ish(iacc)/2.ge.1) then
|
||||
c
|
||||
do 226 i=1,nhops
|
||||
if(lsthop(1,i).eq.isgm(iacc).and.lsthop(2,i).eq.isgm(idon)) then
|
||||
goto 227
|
||||
endif
|
||||
226 continue
|
||||
c
|
||||
c
|
||||
c
|
||||
call qhop_asite(msa,xs,idon,iga(idon),isgm(idon),lseq(isgm(idon)),
|
||||
+ iacc,iga(iacc),isgm(iacc),lseq(isgm(iacc)),xacc,idoff,iaoff)
|
||||
+ iacc,iga(iacc),isgm(iacc),lseq(isgm(iacc)),xacc,idoff,iaoff,
|
||||
+ lsame)
|
||||
c
|
||||
c do not allow hops to the same protonation state if the type of
|
||||
c residue/molecule is identical (e.q. water dissociating into OH and H3O)
|
||||
c
|
||||
if(.not.lsame.or.lseq(isgm(idon)).ne.lseq(isgm(iacc))) then
|
||||
c
|
||||
nlda=nlda+1
|
||||
c write(*,2222) nlda,idon,iga(idon),isgm(idon),lseq(isgm(idon)),
|
||||
c + iacc,iga(iacc),isgm(iacc),lseq(isgm(iacc)),sqrt(rw(jsa))
|
||||
2222 format(9i5,f12.6)
|
||||
c write(*,2222) me,nlda,idon,iga(idon),isgm(idon),lseq(isgm(idon)),
|
||||
c + ish(idon),iacc,iga(iacc),isgm(iacc),lseq(isgm(iacc)),ish(iacc),
|
||||
c + sqrt(rw(jsa))
|
||||
c 2222 format(12i5,f12.6)
|
||||
if(nlda.gt.3*nhop) call md_abort('Increase nhop',nhop)
|
||||
lda(1,nlda)=iga(idon)
|
||||
lda(2,nlda)=idon
|
||||
|
|
@ -619,19 +635,32 @@ c + iacc,iga(iacc),isgm(iacc),lseq(isgm(iacc)),sqrt(rw(jsa))
|
|||
lda(8,nlda)=isgm(iacc)
|
||||
lda(9,nlda)=lseq(isgm(iacc))-ish(iacc)/2
|
||||
lda(10,nlda)=iaoff
|
||||
lda(11,nlda)=idown
|
||||
lda(12,nlda)=0
|
||||
lda(13,nlda)=0
|
||||
lda(14,nlda)=iaown
|
||||
lda(15,nlda)=0
|
||||
lda(16,nlda)=0
|
||||
rda(1,nlda)=xacc(1)
|
||||
rda(2,nlda)=xacc(2)
|
||||
rda(3,nlda)=xacc(3)
|
||||
rda(4,nlda)=xacc(4)
|
||||
rda(5,nlda)=xacc(5)
|
||||
rda(6,nlda)=xacc(6)
|
||||
rda(7,nlda)=xacc(7)
|
||||
rda(8,nlda)=xacc(8)
|
||||
rda(9,nlda)=xacc(9)
|
||||
rda(10,nlda)=xacc(10)
|
||||
rda(11,nlda)=xacc(11)
|
||||
c
|
||||
endif
|
||||
c
|
||||
227 continue
|
||||
c
|
||||
endif
|
||||
endif
|
||||
ctmp
|
||||
endif
|
||||
ctmp
|
||||
c
|
||||
225 continue
|
||||
endif
|
||||
endif
|
||||
|
|
|
|||
|
|
@ -30,7 +30,7 @@ c in/out int : lstptr : list pointer
|
|||
c
|
||||
c dimensions nwm,nwa and nsa need to have been given by a call to cf_initx
|
||||
c
|
||||
c $Id: cf_pairs.F,v 1.14 2003-09-22 22:49:12 d3j191 Exp $
|
||||
c $Id: cf_pairs.F,v 1.15 2006-09-18 15:42:14 d3j191 Exp $
|
||||
implicit none
|
||||
c
|
||||
#include "cf_common.fh"
|
||||
|
|
@ -142,7 +142,8 @@ c
|
|||
+ int_mb(lptr+1+4*nsloc),
|
||||
+ int_mb(i_s2i1),dbl_mb(i_s3r1),dbl_mb(i_s1r1),
|
||||
+ int_mb(i_itrd(2)),mtt(2),int_mb(i_ixcl(2)),mxt(2),
|
||||
+ int_mb(i_lda),dbl_mb(i_rda),lseq,int_mb(i_lsthop))
|
||||
+ int_mb(i_lda),dbl_mb(i_rda),lseq,int_mb(i_lsthop),
|
||||
+ int_mb(i_mprot))
|
||||
endif
|
||||
int_mb(lptr)=npairs
|
||||
int_mb(lptrn+3)=lptr
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
subroutine cf_rdtop(lfntop,filtop,snam)
|
||||
c $Id: cf_rdtop.F,v 1.35 2005-01-08 20:27:06 d3j191 Exp $
|
||||
c $Id: cf_rdtop.F,v 1.36 2006-09-18 15:42:14 d3j191 Exp $
|
||||
implicit none
|
||||
c
|
||||
#include "cf_common.fh"
|
||||
|
|
@ -15,14 +15,14 @@ c
|
|||
real*8 releps
|
||||
character*16 stemp,snam(nsatot)
|
||||
integer natyps,nqtyps
|
||||
real*8 rdata(24),vdata(24,4),qdata(24,2),pdata(4,24)
|
||||
real*8 rdata(24),vdata(24,4),qdata(24,3),pdata(4,24)
|
||||
integer idata(24)
|
||||
character*6 cdata(24)
|
||||
integer i,ia,iq,j,k,l,m,m2,m3,m4,itemp(24),i_i,i_j,l_i,l_j
|
||||
integer naw,nbw,nhw,ndw,now,ntw,nnw
|
||||
integer nhp
|
||||
integer nas,nbt,nhs,nds,nos,nts,nxs
|
||||
integer numbl,numbi,numbd,iqfr,iqto,nsad
|
||||
integer numbl,numbi,numbd,iqfr,iqto,nsad,maxpro
|
||||
character*10 string
|
||||
c
|
||||
real*8 cf_charge
|
||||
|
|
@ -89,16 +89,17 @@ c
|
|||
endif
|
||||
do 4 i=1,nqtyps
|
||||
read(lfntop,1007,end=9997,err=9998)
|
||||
+ ((qdata(k,l),l=1,2),k=1,nparms)
|
||||
1007 format(5x,f12.6,e12.5)
|
||||
+ ((qdata(k,l),l=1,3),k=1,nparms)
|
||||
1007 format(5x,f12.6,e12.5,f12.6)
|
||||
call cf_parq(i,qdata)
|
||||
4 continue
|
||||
if(.not.ma_verify_allocator_stuff()) then
|
||||
call md_abort('ERROR IN MEMORY IN RDTOP CHARGE TYPES',0)
|
||||
endif
|
||||
do 4004 i=1,nseq
|
||||
read(lfntop,4008,end=9997,err=9998) string
|
||||
4008 format(a10)
|
||||
read(lfntop,4008,end=9997,err=9998) string,maxpro
|
||||
4008 format(a10,33x,i5)
|
||||
call cf_parseq(i,maxpro)
|
||||
4004 continue
|
||||
read(lfntop,1008,end=9997,err=9998) naw,nbw,nhw,ndw,now,ntw,nnw
|
||||
1008 format(5i7,2i10)
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
subroutine cf_shakew(xw,yw,iwgm,iwdt,numw)
|
||||
c $Id: cf_shake.F,v 1.26 2003-10-19 03:30:56 d3j191 Exp $
|
||||
c $Id: cf_shake.F,v 1.27 2006-09-18 15:42:15 d3j191 Exp $
|
||||
implicit none
|
||||
c
|
||||
#include "cf_common.fh"
|
||||
|
|
@ -514,7 +514,21 @@ c
|
|||
endif
|
||||
endif
|
||||
endif
|
||||
c
|
||||
c place dummy hydrogens (acceptor sites) onto bound heavy atoms
|
||||
c
|
||||
if(iand(ishp(isa),1).eq.1.and.iand(ishp(jsa),1).ne.1) then
|
||||
xs(isa,1)=xs(jsa,1)
|
||||
xs(isa,2)=xs(jsa,2)
|
||||
xs(isa,3)=xs(jsa,3)
|
||||
endif
|
||||
if(iand(ishp(isa),1).ne.1.and.iand(ishp(jsa),1).eq.1) then
|
||||
xs(jsa,1)=xs(isa,1)
|
||||
xs(jsa,2)=xs(isa,2)
|
||||
xs(jsa,3)=xs(isa,3)
|
||||
endif
|
||||
endif
|
||||
c
|
||||
2 continue
|
||||
c
|
||||
niters=niters+1
|
||||
|
|
|
|||
|
|
@ -1,6 +1,6 @@
|
|||
subroutine md_main()
|
||||
c
|
||||
c $Id: md_main.F,v 1.211 2006-06-05 00:12:29 marat Exp $
|
||||
c $Id: md_main.F,v 1.212 2006-09-18 15:42:03 d3j191 Exp $
|
||||
c
|
||||
implicit none
|
||||
c
|
||||
|
|
@ -97,19 +97,17 @@ c
|
|||
+ epot,epotw,epotsw,epots,volume,dwr,dbl_mb(i_dsr),ekin,etot,
|
||||
+ npolit,dbl_mb(i_gsm),dbl_mb(i_esa),box,dbl_mb(i_xsm))
|
||||
c
|
||||
if(.not.lqmmm)
|
||||
> call prp_print()
|
||||
if(.not.lqmmm) call prp_print()
|
||||
c
|
||||
call rtdb_put(irtdb,'md:energy',mt_dbl,1,epot)
|
||||
c
|
||||
c print energies
|
||||
c
|
||||
if(.not.lqmmm)
|
||||
> call cf_print_energy(lfnout)
|
||||
c
|
||||
if(.not.lqmmm)
|
||||
> call sp_printf(filtop,lfntop,
|
||||
if(.not.lqmmm) then
|
||||
call cf_print_energy(lfnout)
|
||||
call sp_printf(filtop,lfntop,
|
||||
+ int_mb(i_is),dbl_mb(i_xs),dbl_mb(i_fs),npener,dbl_mb(i_esa))
|
||||
endif
|
||||
c
|
||||
if(ifidi.ne.0) then
|
||||
call md_fd(int_mb(i_is),dbl_mb(i_xs),dbl_mb(i_vs),dbl_mb(i_fs),
|
||||
|
|
@ -2481,8 +2479,8 @@ c
|
|||
endif
|
||||
c
|
||||
if(.not.ldone) then
|
||||
if(cf_hopping(stime,isl(1,lssgm),isl(1,lsgan),isl(1,lshop),
|
||||
+ xs,nsaloc,lfnhop)) then
|
||||
if(cf_hopping(lpbc,lpbcs,stime,isl,isl(1,lssgm),isl(1,lsgan),
|
||||
+ isl(1,lsct3),isl(1,lshop),xs,nsaloc,lfnhop)) then
|
||||
if(lpair) then
|
||||
lpair=.false.
|
||||
lforce=.true.
|
||||
|
|
@ -2902,7 +2900,7 @@ c
|
|||
end
|
||||
subroutine md_test
|
||||
c
|
||||
c $Id: md_main.F,v 1.211 2006-06-05 00:12:29 marat Exp $
|
||||
c $Id: md_main.F,v 1.212 2006-09-18 15:42:03 d3j191 Exp $
|
||||
c
|
||||
implicit none
|
||||
c
|
||||
|
|
|
|||
|
|
@ -1,7 +1,7 @@
|
|||
subroutine qhop_asite(msa,xs,idon,idong,resd,istatd,iacc,
|
||||
+ iaccg,resa,istata,xacc,doff,aoff)
|
||||
+ iaccg,resa,istata,xacc,doff,aoff,flada)
|
||||
c
|
||||
c $Id: qhop_asite.F,v 1.3 2003-08-11 19:43:18 d3j191 Exp $
|
||||
c $Id: qhop_asite.F,v 1.4 2006-09-18 15:45:00 d3j191 Exp $
|
||||
c
|
||||
implicit none
|
||||
c
|
||||
|
|
@ -11,11 +11,11 @@ c
|
|||
c declaration of arguments
|
||||
c
|
||||
integer msa,idon,idong,resd,iacc,iaccg,resa,istatd,istata,doff,
|
||||
+ aoff
|
||||
real*8 xacc(5),xs(msa,3)
|
||||
+ aoff,flada
|
||||
real*8 xacc(11),xs(msa,3)
|
||||
c
|
||||
call qhop_asite2(msa,xs,idon,idong,resd,istatd,iacc,
|
||||
+ iaccg,resa,istata,xacc,doff,aoff,
|
||||
+ iaccg,resa,istata,xacc,doff,aoff,flada,
|
||||
+ int_mb(i_ptseq),int_mb(i_ptarat),int_mb(i_iarat),
|
||||
+ dbl_mb(i_racs),dbl_mb(i_racs+mxar),int_mb(i_ptpar),
|
||||
+ dbl_mb(i_par),dbl_mb(i_par+7*maxpar),dbl_mb(i_par+14*maxpar),
|
||||
|
|
@ -25,11 +25,11 @@ c
|
|||
return
|
||||
end
|
||||
subroutine qhop_asite2(msa,xs,idon,idong,resd,istatd,iacc,
|
||||
+ iaccg,resa,istata,xacc,doff,aoff,
|
||||
+ iaccg,resa,istata,xacc,doff,aoff,flada,
|
||||
+ ptseq,ptarat,arat,deq,aneq,ptpar,tunnel,defe12,zpef,
|
||||
+ tdsgl,tstval,e12fxy)
|
||||
c
|
||||
c $Id: qhop_asite.F,v 1.3 2003-08-11 19:43:18 d3j191 Exp $
|
||||
c $Id: qhop_asite.F,v 1.4 2006-09-18 15:45:00 d3j191 Exp $
|
||||
c
|
||||
implicit none
|
||||
c
|
||||
|
|
@ -38,8 +38,8 @@ c
|
|||
c declaration of arguments
|
||||
c
|
||||
integer msa,idon,idong,resd,iacc,iaccg,resa,istatd,istata,doff,
|
||||
+ aoff
|
||||
real*8 xacc(5),xs(msa,3)
|
||||
+ aoff,flada
|
||||
real*8 xacc(11),xs(msa,3)
|
||||
c
|
||||
real*8 tunnel(maxpar,7),defe12(maxpar,7),zpef(maxpar,3),
|
||||
+ tdsgl(maxpar,6),tstval(maxpar,5),e12fxy(maxpar,3)
|
||||
|
|
@ -49,15 +49,15 @@ c
|
|||
c declaration of variables
|
||||
c
|
||||
integer typ,i,idr(2),idhv(2),inh,inhv,numh,nhv(2)
|
||||
real*8 dist,ang
|
||||
real*8 cost,sint,cosp,sinp,th,ph
|
||||
real*8 dist,ang,d2
|
||||
real*8 cost,sint,cosp,sinp,th,ph,alpha
|
||||
real*8 trsl,rds,fct
|
||||
real*8 xh2(3),xacc2(2),xdr(3),xr12(3),xar(3),xdrp(2),norm,r1(3),
|
||||
+ sid1(3),sid2(3),nsid(2),xrh(3),nrh
|
||||
+ sid1(3),sid2(3),nsid(2),xrh(3),nrh,r1pbc(3)
|
||||
real*8 tol,pi
|
||||
c
|
||||
parameter(tol=1E-10)
|
||||
parameter(pi=3.14159)
|
||||
parameter(tol=1E-15)
|
||||
parameter(pi=3.14159265)
|
||||
c
|
||||
typ=arat(ptarat(ptseq(resa,1))+iaccg-
|
||||
+ ptseq(resa,2),2)
|
||||
|
|
@ -72,7 +72,10 @@ c
|
|||
+ ptseq(resa,2))
|
||||
dist=deq(ptarat(ptseq(resa,1))+iaccg-
|
||||
+ ptseq(resa,2))
|
||||
c write (*,98),'angdi ',ang,dist
|
||||
c 98 format(a6,2f10.6)
|
||||
c
|
||||
ang=1.910610
|
||||
idhv(1)=arat(ptarat(ptseq(resd,1))+idong-
|
||||
+ ptseq(resd,2),1)+idon
|
||||
idhv(2)=arat(ptarat(ptseq(resa,1))+iaccg-
|
||||
|
|
@ -102,21 +105,13 @@ c
|
|||
rds=dist*sin(pi-ang)
|
||||
c
|
||||
if(xdr(1).lt.tol.and.xdr(2).lt.tol)then
|
||||
c
|
||||
c no rotation needed
|
||||
c
|
||||
fct=sqrt(xh2(1)**2+xh2(2)**2)
|
||||
if(fct.gt.tol)then
|
||||
xacc(1)=xs(iacc,1)*(rds/fct)
|
||||
xacc(2)=xs(iacc,2)*(rds/fct)
|
||||
else
|
||||
c
|
||||
c H projection on the circle plane coincides with the center of
|
||||
c the circle, all points are equally distant
|
||||
c
|
||||
xacc(1)=rds
|
||||
xacc(2)=0.
|
||||
endif
|
||||
xacc(1)=rds+xs(idhv(2),1)
|
||||
xacc(2)=xs(idhv(2),2)
|
||||
xacc(3)=xs(idhv(2),3)+sign(trsl,xs(idhv(2),3)-xs(idr(1),3))
|
||||
return
|
||||
else
|
||||
ph=atan(-(xdr(2)/(xdr(1))))
|
||||
|
|
@ -196,8 +191,126 @@ c
|
|||
c
|
||||
c H2O: idr(1) and idr(2) are the indexes of the two non dummy H
|
||||
c
|
||||
else
|
||||
c
|
||||
else if(typ.eq.4)then
|
||||
c
|
||||
c if acceptor is water then fix exactly the 2 angles
|
||||
c
|
||||
xdr(1)=xs(idhv(2)+1,1)-xs(idhv(2),1)
|
||||
xdr(2)=xs(idhv(2)+1,2)-xs(idhv(2),2)
|
||||
xdr(3)=xs(idhv(2)+1,3)-xs(idhv(2),3)
|
||||
call cf_pbc(1,xdr,1,xdr,1,0,1,1)
|
||||
c write (*,91),'xdr1 ',xdr
|
||||
|
||||
if(abs(xdr(1)).gt.tol)then
|
||||
if(abs(xdr(2)).gt.tol)then
|
||||
th=atan(-xdr(1)/xdr(2))
|
||||
c write (*,92),'theta ',th,xdr(1),xdr(2)
|
||||
92 format(a6,3f10.6)
|
||||
else
|
||||
th=-sign(1.,xdr(1))*(pi/2)
|
||||
c write *,'thetanotol ',xdr(1)-tol,xdr(2)-tol
|
||||
endif
|
||||
else
|
||||
th=0
|
||||
endif
|
||||
xdrp(1)=-xdr(1)*sin(th)+xdr(2)*cos(th)
|
||||
|
||||
if(abs(xdrp(1)).gt.tol)then
|
||||
if(abs(xdr(3)).gt.tol)then
|
||||
ph=atan(-xdrp(1)/xdr(3))
|
||||
c write (*,92),'phi ',ph,xdrp(1),xdr(3)
|
||||
else
|
||||
ph=-sign(1.,xdrp(1))
|
||||
c write *,'phinotol ',xdrp(1)-tol,xdr(3)-tol
|
||||
endif
|
||||
else
|
||||
ph=0
|
||||
endif
|
||||
xdrp(2)=-xdrp(1)*sin(ph)+xdr(3)*cos(ph)
|
||||
if(xdrp(2).le.0)ph=ph+pi
|
||||
|
||||
xdr(1)=xs(idhv(2)+2,1)-xs(idhv(2),1)
|
||||
xdr(2)=xs(idhv(2)+2,2)-xs(idhv(2),2)
|
||||
xdr(3)=xs(idhv(2)+2,3)-xs(idhv(2),3)
|
||||
call cf_pbc(1,xdr,1,xdr,1,0,1,1)
|
||||
|
||||
call rot(th,xdr(1),xdr(2),xdr(3),3)
|
||||
call rot(ph,xdr(1),xdr(2),xdr(3),1)
|
||||
c write (*,91),'xdr2 ',xdr
|
||||
if(abs(xdr(2)).gt.tol)then
|
||||
if(abs(xdr(1)).gt.tol)then
|
||||
alpha=atan(xdr(2)/xdr(1))
|
||||
c write (*,92),'alpha ',alpha,xdr(2),xdr(1)
|
||||
else
|
||||
alpha=sign(1.,xdr(2))
|
||||
endif
|
||||
else
|
||||
alpha=0
|
||||
endif
|
||||
xdrp(1)=xdr(1)*cos(alpha)+xdr(2)*sin(alpha)
|
||||
if(xdrp(1).le.0)then
|
||||
if(alpha.le.0)then
|
||||
alpha=alpha+pi
|
||||
else
|
||||
alpha=alpha-pi
|
||||
endif
|
||||
endif
|
||||
c
|
||||
c Found the angles, perform rotation
|
||||
c
|
||||
xdr(1)=xs(idhv(1),1)-xs(idhv(2),1)
|
||||
xdr(2)=xs(idhv(1),2)-xs(idhv(2),2)
|
||||
xdr(3)=xs(idhv(1),3)-xs(idhv(2),3)
|
||||
call cf_pbc(1,xdr,1,xdr,1,0,1,1)
|
||||
|
||||
call rot(th,xdr(1),xdr(2),xdr(3),3)
|
||||
call rot(ph,xdr(1),xdr(2),xdr(3),1)
|
||||
call rot(alpha,xdr(1),xdr(2),xdr(3),3)
|
||||
xr12(1)=xdr(1)
|
||||
xr12(2)=xdr(2)
|
||||
xr12(3)=xdr(3)
|
||||
xdr(1)=xs(idhv(2)+2,1)-xs(idhv(2),1)
|
||||
xdr(2)=xs(idhv(2)+2,2)-xs(idhv(2),2)
|
||||
xdr(3)=xs(idhv(2)+2,3)-xs(idhv(2),3)
|
||||
call cf_pbc(1,xdr,1,xdr,1,0,1,1)
|
||||
|
||||
call rot(th,xdr(1),xdr(2),xdr(3),3)
|
||||
call rot(ph,xdr(1),xdr(2),xdr(3),1)
|
||||
call rot(alpha,xdr(1),xdr(2),xdr(3),3)
|
||||
c write (*,91),'xdr3 ',xdr
|
||||
d2=dist**2
|
||||
xacc(3)=dist*cos(ang)
|
||||
xacc(1)=(1/xdr(1))*xacc(3)*(0.1-xdr(3))
|
||||
xacc(2)=sqrt(d2-xacc(3)**2-xacc(1)**2)
|
||||
c write (*,93),'xacc1 ',xacc(1),xacc(2),xacc(3),d2,cos(ang)
|
||||
c write(*,94),'norm ',xacc(1)**2+xacc(2)**2+xacc(3)**2
|
||||
d2=abs(xr12(2)-xacc(2))-abs(xr12(2)+xacc(2))
|
||||
if(d2.ge.0)xacc(2)=-xacc(2)
|
||||
xdr(1)=xacc(1)
|
||||
xdr(2)=xacc(2)
|
||||
xdr(3)=xacc(3)
|
||||
call rot(-alpha,xdr(1),xdr(2),xdr(3),3)
|
||||
call rot(-ph,xdr(1),xdr(2),xdr(3),1)
|
||||
call rot(-th,xdr(1),xdr(2),xdr(3),3)
|
||||
c write (*,94),'norm ',xacc(1)**2+xacc(2)**2+xacc(3)**2
|
||||
xacc(1)=xdr(1)+xs(idhv(2),1)
|
||||
xacc(2)=xdr(2)+xs(idhv(2),2)
|
||||
xacc(3)=xdr(3)+xs(idhv(2),3)
|
||||
c write (*,93),'xacc2 ',xacc(1),xacc(2),xacc(3),d2,cos(ang)
|
||||
91 format(a5,3f10.6)
|
||||
93 format(a6,5f10.6)
|
||||
94 format(a5,f10.6)
|
||||
|
||||
else if(typ.eq.5)then
|
||||
|
||||
xar(1)=xs(idon,1)-xs(idhv(2),1)
|
||||
xar(2)=xs(idon,2)-xs(idhv(2),2)
|
||||
xar(3)=xs(idon,3)-xs(idhv(2),3)
|
||||
call cf_pbc(1,xar,1,xar,1,0,1,1)
|
||||
norm=sqrt(xar(1)**2+xar(2)**2+xar(3)**2)
|
||||
xacc(1)=xs(idhv(2),1)+dist*(xar(1)/norm)
|
||||
xacc(2)=xs(idhv(2),2)+dist*(xar(2)/norm)
|
||||
xacc(3)=xs(idhv(2),3)+dist*(xar(3)/norm)
|
||||
endif
|
||||
c
|
||||
c Calculating angle and donor-acceptor dist
|
||||
|
|
@ -205,8 +318,37 @@ c
|
|||
r1(1)=xs(idhv(1),1)-xs(idhv(2),1)
|
||||
r1(2)=xs(idhv(1),2)-xs(idhv(2),2)
|
||||
r1(3)=xs(idhv(1),3)-xs(idhv(2),3)
|
||||
r1pbc(1)=r1(1)
|
||||
r1pbc(2)=r1(2)
|
||||
r1pbc(3)=r1(3)
|
||||
call cf_pbc(1,r1,1,r1,1,0,1,1)
|
||||
xacc(4)=sqrt(r1(1)**2+r1(2)**2+r1(3)**2)
|
||||
c
|
||||
c calculating midpoint of donor-acceptor
|
||||
c
|
||||
if(abs(r1pbc(1)-r1(1)).lt.0.001) then
|
||||
xacc(6)=(xs(idhv(1),1)+xs(idhv(2),1))/2
|
||||
xacc(9)=xacc(6)
|
||||
else
|
||||
xacc(6)=(xs(idhv(1),1)+(r1pbc(1)-r1(1))+xs(idhv(2),1))/2
|
||||
xacc(9)=xacc(6)-(r1pbc(1)-r1(1))
|
||||
endif
|
||||
c
|
||||
if(abs(r1pbc(2)-r1(2)).lt.0.001) then
|
||||
xacc(7)=(xs(idhv(1),2)+xs(idhv(2),2))/2
|
||||
xacc(10)=xacc(7)
|
||||
else
|
||||
xacc(7)=(xs(idhv(1),2)+(r1pbc(2)-r1(2))+xs(idhv(2),2))/2
|
||||
xacc(10)=xacc(7)-(r1pbc(2)-r1(2))
|
||||
endif
|
||||
c
|
||||
if(abs(r1pbc(3)-r1(3)).lt.0.001) then
|
||||
xacc(8)=(xs(idhv(1),3)+xs(idhv(2),3))/2
|
||||
xacc(11)=xacc(8)
|
||||
else
|
||||
xacc(8)=(xs(idhv(1),3)+(r1pbc(3)-r1(3))+xs(idhv(2),3))/2
|
||||
xacc(11)=xacc(8)-(r1pbc(3)-r1(3))
|
||||
endif
|
||||
c
|
||||
sid1(1)=xs(idhv(1),1)-xs(idon,1)
|
||||
sid1(2)=xs(idhv(1),2)-xs(idon,2)
|
||||
|
|
@ -224,6 +366,15 @@ c
|
|||
+ (1/(nsid(1)*nsid(2)))
|
||||
xacc(5)=acos(xacc(5))
|
||||
c
|
||||
c
|
||||
c debug: print shift of proton
|
||||
c
|
||||
r1(1)=xs(idon,1)-xacc(1)
|
||||
r1(2)=xs(idon,2)-xacc(2)
|
||||
r1(3)=xs(idon,3)-xacc(3)
|
||||
call cf_pbc(1,r1,1,r1,1,0,1,1)
|
||||
c dshop=sqrt(r1(1)**2+r1(2)**2+r1(3)**2)
|
||||
c print*,"mov ",resd,idong,resa,norm,xacc(4),xacc(5)
|
||||
c Shift for permutation of donor H
|
||||
c
|
||||
numh=arat(ptarat(ptseq(resd,1))+idong-
|
||||
|
|
@ -232,8 +383,15 @@ c
|
|||
+ ptseq(resd,2),1)
|
||||
inhv=arat(ptarat(ptseq(resd,1))+idong-
|
||||
+ ptseq(resd,2),3)
|
||||
doff=numh-inh-iand(((istatd-1)/(2**(inhv-1))+1),1)
|
||||
if(typ.eq.3)doff=numh-inh-(istatd-1)
|
||||
if(numh.eq.1)then
|
||||
doff=0
|
||||
else if(numh.eq.2)then
|
||||
doff=numh-inh
|
||||
c doff=numh-inh-iand(((istatd-1)/(2**(inhv-1))+1),1)
|
||||
else
|
||||
doff=numh-inh-(istatd-1)
|
||||
endif
|
||||
c print*,'don ',resd,inh,doff,istatd
|
||||
c
|
||||
c Shift for permutation of acceptor dummy site
|
||||
c Currently if nhv>1,only one dummy is allowed
|
||||
|
|
@ -246,11 +404,39 @@ c
|
|||
+ ptseq(resa,2),3)
|
||||
aoff=inh-numh-1+iand(((istata-1)/(2**(inhv-1))+1),1)
|
||||
aoff=0
|
||||
if(typ.eq.3)aoff=numh-inh-(istata-2)
|
||||
c print*,numh,inh,istata,aoff
|
||||
c if(numh>2)aoff=numh-inh-(istata-2)
|
||||
c print*,'acc ',resa,inh,doff,istata
|
||||
c
|
||||
c print 1021,typ,doff,aoff,xacc,numh,inh,istatd
|
||||
1021 format(3i5,5f8.3,3i5)
|
||||
c 1021 format(3i5,5f8.3,3i5)
|
||||
c
|
||||
c flag to see if res are same type
|
||||
c
|
||||
flada=0
|
||||
if(ptseq(resa,1).eq.ptseq(resd,1))flada=1
|
||||
return
|
||||
end
|
||||
|
||||
subroutine rot(angle,x,y,z,axis)
|
||||
|
||||
implicit none
|
||||
|
||||
real*8 angle,x,y,z,x1,y1,z1
|
||||
integer axis
|
||||
|
||||
x1=x
|
||||
y1=y
|
||||
z1=z
|
||||
|
||||
if(axis.eq.3)then
|
||||
x=x1*cos(angle)+y1*sin(angle)
|
||||
y=-x1*sin(angle)+y1*cos(angle)
|
||||
else if(axis.eq.1) then
|
||||
y=y1*cos(angle)+z1*sin(angle)
|
||||
z=-y1*sin(angle)+z1*cos(angle)
|
||||
else
|
||||
|
||||
endif
|
||||
end
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -1,4 +1,4 @@
|
|||
subroutine qhop_prob(gd,rsd,prd,ga,rsa,pra,uda,dang,prob)
|
||||
subroutine qhop_prob(gd,rsd,prd,ga,rsa,pra,uda,dang,prob,e120)
|
||||
c
|
||||
implicit none
|
||||
c
|
||||
|
|
@ -7,7 +7,7 @@ c
|
|||
c
|
||||
c declaration of arguments
|
||||
c
|
||||
real*8 uda(4),dang(2),prob
|
||||
real*8 uda(4),dang(2),prob,e120
|
||||
integer gd,rsd,prd,ga,rsa,pra
|
||||
c
|
||||
call qhop_prob2(gd,rsd,prd,ga,rsa,pra,uda,dang,prob,
|
||||
|
|
@ -15,15 +15,17 @@ c
|
|||
+ dbl_mb(i_racs),dbl_mb(i_racs+mxar),int_mb(i_ptpar),
|
||||
+ dbl_mb(i_par),dbl_mb(i_par+7*maxpar),dbl_mb(i_par+14*maxpar),
|
||||
+ dbl_mb(i_par+17*maxpar),dbl_mb(i_par+23*maxpar),
|
||||
+ dbl_mb(i_par+28*maxpar))
|
||||
+ dbl_mb(i_par+28*maxpar),e120)
|
||||
c
|
||||
e120=e120*4.184
|
||||
c
|
||||
return
|
||||
end
|
||||
subroutine qhop_prob2(gd,rsd,prd,ga,rsa,pra,uda,dang,prob,
|
||||
+ ptseq,ptarat,arat,deq,aneq,ptpar,tunnel,defe12,zpef,
|
||||
+ tdsgl,tstval,e12fxy)
|
||||
+ tdsgl,tstval,e12fxy,e120)
|
||||
c
|
||||
c $Id: qhop_prob.F,v 1.6 2003-09-22 15:20:10 d3j191 Exp $
|
||||
c $Id: qhop_prob.F,v 1.7 2006-09-18 15:45:00 d3j191 Exp $
|
||||
c
|
||||
c
|
||||
c The hopping probability prob is calculated given the difference in
|
||||
|
|
@ -46,7 +48,7 @@ c
|
|||
c
|
||||
c declaration of arguments
|
||||
c
|
||||
real*8 uda(4),dang(2),prob
|
||||
real*8 uda(4),dang(2),prob,e120
|
||||
integer gd,rsd,prd,ga,rsa,pra
|
||||
c
|
||||
real*8 tunnel(maxpar,7),defe12(maxpar,7),zpef(maxpar,3),
|
||||
|
|
@ -100,9 +102,8 @@ c 1026 format('prob ',9i5)
|
|||
c
|
||||
c Calculating e12
|
||||
c
|
||||
e12=(uda(3)+uda(4)-uda(1)-uda(2))/4.184
|
||||
e12=e12+e12fxy(ipar,1)+e12fxy(ipar,2)*dist+
|
||||
+ e12fxy(ipar,3)*(dist**2)
|
||||
e120=e12fxy(ipar,1)+e12fxy(ipar,2)*dist+e12fxy(ipar,3)*(dist**2)
|
||||
e12=e12+e120+(uda(3)+uda(4)-uda(1)-uda(2))/4.184
|
||||
c
|
||||
c lhv(1)=ld+arat(ptarat(ptseq(iseq(1),1)+gd-ptseq(iseq(1),2),1)
|
||||
c lhv(2)=la+arat(ptarat(ptseq(iseq(2),1)+ga-ptseq(iseq(2),2),1)
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
#
|
||||
# $Id: GNUmakefile,v 1.3 2006-08-25 20:24:45 d3p307 Exp $
|
||||
# $Id: GNUmakefile,v 1.4 2006-09-18 15:44:09 d3j191 Exp $
|
||||
#
|
||||
#
|
||||
OBJ_OPTIMIZE=sp_init.o sp_main.o
|
||||
|
|
@ -16,15 +16,4 @@ ifeq ($(TARGET),LINUX)
|
|||
DEFINES += -DNEED_LOC
|
||||
endif
|
||||
|
||||
ifeq ($(FC),ifc)
|
||||
FOPTIMIZE += -Ob0
|
||||
endif
|
||||
|
||||
ifeq ($(FC),ifort)
|
||||
ifeq ($(_CPU),x86_64)
|
||||
FOPTIMIZE += -O0
|
||||
endif
|
||||
endif
|
||||
|
||||
|
||||
include ../config/makelib.h
|
||||
|
|
|
|||
|
|
@ -1,7 +1,7 @@
|
|||
subroutine sp_setup(node,iwl,xw,xwcr,vw,gw,numwm,
|
||||
+ isl,xs,xscr,vs,gs,numsa,lpack)
|
||||
c
|
||||
c $Id: sp_main.F,v 1.97 2006-03-07 19:44:28 d3j191 Exp $
|
||||
c $Id: sp_main.F,v 1.98 2006-09-18 15:44:09 d3j191 Exp $
|
||||
c
|
||||
implicit none
|
||||
c
|
||||
|
|
@ -2866,6 +2866,65 @@ c
|
|||
+ int_mb(i_bb),int_mb(i_ipl),int_mb(i_jpl),lself,.false.)
|
||||
c
|
||||
lpbcs=lpbc9
|
||||
c
|
||||
return
|
||||
end
|
||||
subroutine sp_getda(jnode,jsfr,jsto,isl,isga,isq3,xs,ndxda)
|
||||
c
|
||||
implicit none
|
||||
c
|
||||
#include "sp_common.fh"
|
||||
#include "mafdecls.fh"
|
||||
c
|
||||
integer jnode,jsfr,jsto,ndxda
|
||||
integer isga(msa),isq3(msa),isl(msa,mis2)
|
||||
real*8 xs(msa,3)
|
||||
c
|
||||
integer il,ih,jl,jh
|
||||
integer l_islt,i_islt
|
||||
c
|
||||
if(.not.ma_push_get(mt_int,msa*mis2,'islt',l_islt,i_islt))
|
||||
+ call md_abort('Failed to allocate memory for islt',0)
|
||||
call sp_gtda(jnode,jsfr,jsto,int_mb(i_islt),int_mb(i_pack),
|
||||
+ isga,isq3,xs,ndxda)
|
||||
if(.not.ma_pop_stack(l_islt))
|
||||
+ call md_abort('Failed to deallocate memory for islt',0)
|
||||
c
|
||||
return
|
||||
end
|
||||
subroutine sp_gtda(jnode,jsfr,jsto,isl,islp,isga,isq3,xs,ndxda)
|
||||
c
|
||||
implicit none
|
||||
c
|
||||
#include "sp_common.fh"
|
||||
#include "global.fh"
|
||||
c
|
||||
integer jnode,jsfr,jsto,ndxda
|
||||
integer isl(msa,mis2),islp(msa,npack)
|
||||
integer isga(msa),isq3(msa)
|
||||
real*8 xs(msa,3)
|
||||
c
|
||||
integer i,il,ih,jl,jh
|
||||
c
|
||||
if(npack.eq.0) then
|
||||
call ga_distribution(ga_is,jnode,il,ih,jl,jh)
|
||||
call ga_get(ga_is,il+jsfr-1,il+jsto-1,jl,jl,
|
||||
+ isga(ndxda+1),msa)
|
||||
call ga_get(ga_is,il+jsfr-1,il+jsto-1,jl+lsct3-1,jl+lsct3-1,
|
||||
+ isq3(ndxda+1),msa)
|
||||
else
|
||||
call ga_distribution(ga_is,jnode,il,ih,jl,jh)
|
||||
call ga_get(ga_is,il+jsfr-1,il+jsto-1,jl,jl+npack-1,
|
||||
+ islp(nsaloc+1,1),msa)
|
||||
call sp_unpack(jsto-jsfr+1,isl(nsaloc+1,1),islp(nsaloc+1,1))
|
||||
do 1 i=1,jsto-jsfr
|
||||
isga(ndxda+i)=isl(nsaloc+i,lsgan)
|
||||
isq3(ndxda+i)=isl(nsaloc+i,lsct3)
|
||||
1 continue
|
||||
endif
|
||||
call ga_distribution(ga_s,jnode,il,ih,jl,jh)
|
||||
call ga_get(ga_s,il+jsfr-1,il+jsto-1,jl,jl+2,
|
||||
+ xs(ndxda+1,1),msa)
|
||||
c
|
||||
return
|
||||
end
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue