- module energy
+ module energy
!-----------------------------------------------------------------------------
use io_units
use names
! amino-acid residue.
! common /precomp1/
real(kind=8),dimension(:,:),allocatable :: mu,muder,Ub2,Ub2der,&
- Ctobr,Ctobrder,Dtobr2,Dtobr2der !(2,maxres)
+ Ctobr,Ctobrder,Dtobr2,Dtobr2der,gUb2 !(2,maxres)
real(kind=8),dimension(:,:,:),allocatable :: EUg,EUgder,CUg,&
CUgder,DUg,Dugder,DtUg2,DtUg2der !(2,2,maxres)
! This common block contains vectors and matrices dependent on two
real(kind=8),dimension(:,:,:,:),allocatable :: Ug2DtEUgder,&
DtUg2EUgder !(2,2,2,maxres)
! common /rotat_old/
+ real(kind=8),dimension(4) :: gmuij,gmuij1,gmuij2,gmuji1,gmuji2
real(kind=8),dimension(:),allocatable :: costab,sintab,&
costab2,sintab2 !(maxres)
! This common block contains dipole-interaction matrices and their
ebe_nucl,esbloc,etors_nucl,etors_d_nucl,ecorr_nucl,&
ecorr3_nucl
! energies for ions
- real(kind=8) :: ecation_prot,ecationcation
+ real(kind=8) :: ecation_prot,ecationcation,ecations_prot_amber
! energies for protein nucleic acid interaction
real(kind=8) :: escbase,epepbase,escpho,epeppho
! print *,"Processor",myrank," BROADCAST iorder"
! FG master sets up the WEIGHTS_ array which will be broadcast to the
! FG slaves as WEIGHTS array.
- ! weights_(1)=wsc
+ weights_(1)=wsc
weights_(2)=wscp
weights_(3)=welec
weights_(4)=wcorr
weights_(41)=wcatcat
weights_(42)=wcatprot
weights_(46)=wscbase
- weights_(47)=wscpho
- weights_(48)=wpeppho
+ weights_(47)=wpepbase
+ weights_(48)=wscpho
+ weights_(49)=wpeppho
! wcatcat= weights(41)
! wcatprot=weights(42)
wcatcat= weights(41)
wcatprot=weights(42)
wscbase=weights(46)
- wscpho=weights(47)
- wpeppho=weights(48)
+ wpepbase=weights(47)
+ wscpho=weights(48)
+ wpeppho=weights(49)
+! welpsb=weights(28)*fact(1)
+!
+! wcorr_nucl= weights(37)*fact(1)
+! wcorr3_nucl=weights(38)*fact(2)
+! wtor_nucl= weights(35)*fact(1)
+! wtor_d_nucl=weights(36)*fact(2)
+
endif
time_Bcast=time_Bcast+MPI_Wtime()-time00
time_Bcastw=time_Bcastw+MPI_Wtime()-time00
.or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.d0 &
.or. wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0 ) then
#endif
- write(iout,*),"just befor eelec call"
+! print *,"just befor eelec call"
call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
-! write (iout,*) "ELEC calc"
+! print *, "ELEC calc"
else
ees=0.0d0
evdw1=0.0d0
! Calculate the virtual-bond-angle energy.
! write(iout,*) "in etotal afer edis",ipot
- if (wang.gt.0.0d0) then
- call ebend(ebe,ethetacnstr)
+! if (wang.gt.0.0d0) then
+! call ebend(ebe,ethetacnstr)
+! else
+! ebe=0
+! ethetacnstr=0
+! endif
+ if (wang.gt.0d0) then
+ if (tor_mode.eq.0) then
+ call ebend(ebe)
+ else
+!C ebend kcc is Kubo cumulant clustered rigorous attemp to derive the
+!C energy function
+ call ebend_kcc(ebe)
+ endif
else
- ebe=0
- ethetacnstr=0
+ ebe=0.0d0
endif
+ ethetacnstr=0.0d0
+ if (with_theta_constr) call etheta_constr(ethetacnstr)
+
! write(iout,*) "in etotal afer ebe",ipot
! print *,"Processor",myrank," computed UB"
! Calculate the virtual-bond torsional energy.
!
!d print *,'nterm=',nterm
- if (wtor.gt.0) then
- call etor(etors,edihcnstr)
+! if (wtor.gt.0) then
+! call etor(etors,edihcnstr)
+! else
+! etors=0
+! edihcnstr=0
+! endif
+ if (wtor.gt.0.0d0) then
+ if (tor_mode.eq.0) then
+ call etor(etors)
+ else
+!C etor kcc is Kubo cumulant clustered rigorous attemp to derive the
+!C energy function
+ call etor_kcc(etors)
+ endif
else
- etors=0
- edihcnstr=0
+ etors=0.0d0
endif
+ edihcnstr=0.0d0
+ if (ndih_constr.gt.0) call etor_constr(edihcnstr)
+!c print *,"Processor",myrank," computed Utor"
+
! print *,"Processor",myrank," computed Utor"
!
call AFMforce(Eafmforce)
else if (selfguide.gt.0) then
call AFMvel(Eafmforce)
+ else
+ Eafmforce=0.0d0
endif
endif
if (tubemode.eq.1) then
eespp=0.0d0
endif
! write(iout,*) ecorr_nucl,"ecorr_nucl",nres_molec(2)
+! print *,"before ecatcat",wcatcat
if (nfgtasks.gt.1) then
if (fg_rank.eq.0) then
call ecatcat(ecationcation)
call ecatcat(ecationcation)
endif
call ecat_prot(ecation_prot)
+ call ecats_prot_amber(ecations_prot_amber)
if (nres_molec(2).gt.0) then
call eprot_sc_base(escbase)
call epep_sc_base(epepbase)
epeppho=0.0
endif
! call ecatcat(ecationcation)
-! print *,"after ebend", ebe_nucl
+! print *,"after ebend", wtor_nucl
#ifdef TIMING
time_enecalc=time_enecalc+MPI_Wtime()-time00
#endif
! Here are the energies showed per procesor if the are more processors
! per molecule then we sum it up in sum_energy subroutine
! print *," Processor",myrank," calls SUM_ENERGY"
- energia(41)=ecation_prot
- energia(42)=ecationcation
+ energia(42)=ecation_prot
+ energia(41)=ecationcation
energia(46)=escbase
energia(47)=epepbase
energia(48)=escpho
energia(49)=epeppho
+ energia(50)=ecations_prot_amber
call sum_energy(energia,.true.)
if (dyn_ss) call dyn_set_nss
! print *," Processor",myrank," left SUM_ENERGY"
real(kind=8) :: evdwpp,eespp,evdwpsb,eelpsb,evdwsb,eelsb,estr_nucl,&
ebe_nucl,esbloc,etors_nucl,etors_d_nucl,ecorr_nucl,&
ecorr3_nucl
- real(kind=8) :: ecation_prot,ecationcation
+ real(kind=8) :: ecation_prot,ecationcation,ecations_prot_amber
real(kind=8) :: escbase,epepbase,escpho,epeppho
integer :: i
#ifdef MPI
etors_d_nucl=energia(36)
ecorr_nucl=energia(37)
ecorr3_nucl=energia(38)
- ecation_prot=energia(41)
- ecationcation=energia(42)
+ ecation_prot=energia(42)
+ ecationcation=energia(41)
escbase=energia(46)
epepbase=energia(47)
escpho=energia(48)
epeppho=energia(49)
+ ecations_prot_amber=energia(50)
+
! energia(41)=ecation_prot
! energia(42)=ecationcation
+wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
+wcatprot*ecation_prot+wcatcat*ecationcation+wscbase*escbase&
- +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho
+ +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+ecations_prot_amber
#else
etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1) &
+wang*ebe+wtor*etors+wscloc*escloc &
+wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
+wcatprot*ecation_prot+wcatcat*ecationcation+wscbase*escbase&
- +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho
+ +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+ecations_prot_amber
#endif
energia(0)=etot
! detecting NaNQ
wtor=weights(13)*fact(1)
wtor_d=weights(14)*fact(2)
wsccor=weights(21)*fact(1)
-
+ welpsb=weights(28)*fact(1)
+ wcorr_nucl= weights(37)*fact(1)
+ wcorr3_nucl=weights(38)*fact(2)
+ wtor_nucl= weights(35)*fact(1)
+ wtor_d_nucl=weights(36)*fact(2)
+ wpepbase=weights(47)*fact(1)
return
end subroutine rescale_weights
!-----------------------------------------------------------------------------
real(kind=8) :: evdwpp,eespp,evdwpsb,eelpsb,evdwsb,eelsb,estr_nucl,&
ebe_nucl,esbloc,etors_nucl,etors_d_nucl,ecorr_nucl,&
ecorr3_nucl
- real(kind=8) :: ecation_prot,ecationcation
+ real(kind=8) :: ecation_prot,ecationcation,ecations_prot_amber
real(kind=8) :: escbase,epepbase,escpho,epeppho
etot=energia(0)
etors_d_nucl=energia(36)
ecorr_nucl=energia(37)
ecorr3_nucl=energia(38)
- ecation_prot=energia(41)
- ecationcation=energia(42)
+ ecation_prot=energia(42)
+ ecationcation=energia(41)
escbase=energia(46)
epepbase=energia(47)
escpho=energia(48)
epeppho=energia(49)
+ ecations_prot_amber=energia(50)
#ifdef SPLITELE
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,&
estr,wbond,ebe,wang,&
etors_d_nucl,wtor_d_nucl,ecorr_nucl,wcorr_nucl,&
ecorr3_nucl,wcorr3_nucl,ecation_prot,wcatprot,ecationcation,wcatcat, &
escbase,wscbase,epepbase,wpepbase,escpho,wscpho,epeppho,wpeppho,&
- etot
+ ecations_prot_amber,etot
10 format (/'Virtual-chain energies:'// &
'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ &
'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ &
ecorr,wcorr,&
ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,&
eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,edihcnstr,&
- ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforc, &
+ ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforce, &
etube,wtube, &
estr_nucl,wbond_nucl, ebe_nucl,wang_nucl,&
- evdwpp,wvdwpp_nucl,eespp,welpp,evdwpsb,wvdwpsb,eelpsb,welpsb&
- evdwsb,wvdwsb,eelsb,welsb,esbloc,wsbloc,etors_nucl,wtor_nucl&
+ evdwpp,wvdwpp_nucl,eespp,welpp,evdwpsb,wvdwpsb,eelpsb,welpsb,&
+ evdwsb,wvdwsb,eelsb,welsb,esbloc,wsbloc,etors_nucl,wtor_nucl,&
etors_d_nucl,wtor_d_nucl,ecorr_nucl,wcorr_nucl,&
ecorr3_nucl,wcorr3_nucl,ecation_prot,wcatprot,ecationcation,wcatcat, &
escbase,wscbase,epepbase,wpepbase,escpho,wscpho,epeppho,wpeppho,&
- etot
+ ecations_prot_amber,etot
10 format (/'Virtual-chain energies:'// &
'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ &
'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ &
! write(iout,*)"c ", c(1,:), c(2,:), c(3,:)
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
rij=dsqrt(rrij)
- sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj)))
- sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj)))
+ sss_ele_cut=sscale_ele(1.0d0/(rij))
+ sss_ele_grad=sscagrad_ele(1.0d0/(rij))
! print *,sss_ele_cut,sss_ele_grad,&
! 1.0d0/(rij),r_cut_ele,rlamb_ele
if (sss_ele_cut.le.0.0) cycle
fac=rij*fac
! print *,'before fac',fac,rij,evdwij
fac=fac+evdwij*sss_ele_grad/sss_ele_cut&
- /sigma(itypi,itypj)*rij
+ *rij
! print *,'grad part scale',fac, &
! evdwij*sss_ele_grad/sss_ele_cut &
! /sigma(itypi,itypj)*rij
! include 'COMMON.FFIELD'
real(kind=8) :: auxvec(2),auxmat(2,2)
integer :: i,iti1,iti,k,l
- real(kind=8) :: sin1,cos1,sin2,cos2,dwacos2,dwasin2
+ real(kind=8) :: sin1,cos1,sin2,cos2,dwacos2,dwasin2,cost1,sint1,&
+ sint1sq,sint1cub,sint1cost1,b1k,b2k,aux
! print *,"in set matrices"
!
! Compute the virtual-bond-torsional-angle dependent quantities needed
! to calculate the el-loc multibody terms of various order.
!
!AL el mu=0.0d0
+
+#ifdef PARMAT
+ do i=ivec_start+2,ivec_end+2
+#else
+ do i=3,nres+1
+#endif
+ if (i.gt. nnt+2 .and. i.lt.nct+2) then
+ if (itype(i-2,1).eq.0) then
+ iti = nloctyp
+ else
+ iti = itype2loc(itype(i-2,1))
+ endif
+ else
+ iti=nloctyp
+ endif
+!c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
+ if (i.gt. nnt+1 .and. i.lt.nct+1) then
+ iti1 = itype2loc(itype(i-1,1))
+ else
+ iti1=nloctyp
+ endif
+! print *,i,itype(i-2,1),iti
+#ifdef NEWCORR
+ cost1=dcos(theta(i-1))
+ sint1=dsin(theta(i-1))
+ sint1sq=sint1*sint1
+ sint1cub=sint1sq*sint1
+ sint1cost1=2*sint1*cost1
+! print *,"cost1",cost1,theta(i-1)
+!c write (iout,*) "bnew1",i,iti
+!c write (iout,*) (bnew1(k,1,iti),k=1,3)
+!c write (iout,*) (bnew1(k,2,iti),k=1,3)
+!c write (iout,*) "bnew2",i,iti
+!c write (iout,*) (bnew2(k,1,iti),k=1,3)
+!c write (iout,*) (bnew2(k,2,iti),k=1,3)
+ k=1
+! print *,bnew1(1,k,iti),"bnew1"
+ do k=1,2
+ b1k=bnew1(1,k,iti)+(bnew1(2,k,iti)+bnew1(3,k,iti)*cost1)*cost1
+! print *,b1k
+! write(*,*) shape(b1)
+! if(.not.allocated(b1)) print *, "WTF?"
+ b1(k,i-2)=sint1*b1k
+!
+! print *,b1(k,i-2)
+
+ gtb1(k,i-2)=cost1*b1k-sint1sq*&
+ (bnew1(2,k,iti)+2*bnew1(3,k,iti)*cost1)
+! print *,gtb1(k,i-2)
+
+ b2k=bnew2(1,k,iti)+(bnew2(2,k,iti)+bnew2(3,k,iti)*cost1)*cost1
+ b2(k,i-2)=sint1*b2k
+! print *,b2(k,i-2)
+
+ gtb2(k,i-2)=cost1*b2k-sint1sq*&
+ (bnew2(2,k,iti)+2*bnew2(3,k,iti)*cost1)
+! print *,gtb2(k,i-2)
+
+ enddo
+! print *,b1k,b2k
+ do k=1,2
+ aux=ccnew(1,k,iti)+(ccnew(2,k,iti)+ccnew(3,k,iti)*cost1)*cost1
+ cc(1,k,i-2)=sint1sq*aux
+ gtcc(1,k,i-2)=sint1cost1*aux-sint1cub*&
+ (ccnew(2,k,iti)+2*ccnew(3,k,iti)*cost1)
+ aux=ddnew(1,k,iti)+(ddnew(2,k,iti)+ddnew(3,k,iti)*cost1)*cost1
+ dd(1,k,i-2)=sint1sq*aux
+ gtdd(1,k,i-2)=sint1cost1*aux-sint1cub*&
+ (ddnew(2,k,iti)+2*ddnew(3,k,iti)*cost1)
+ enddo
+! print *,"after cc"
+ cc(2,1,i-2)=cc(1,2,i-2)
+ cc(2,2,i-2)=-cc(1,1,i-2)
+ gtcc(2,1,i-2)=gtcc(1,2,i-2)
+ gtcc(2,2,i-2)=-gtcc(1,1,i-2)
+ dd(2,1,i-2)=dd(1,2,i-2)
+ dd(2,2,i-2)=-dd(1,1,i-2)
+ gtdd(2,1,i-2)=gtdd(1,2,i-2)
+ gtdd(2,2,i-2)=-gtdd(1,1,i-2)
+! print *,"after dd"
+
+ do k=1,2
+ do l=1,2
+ aux=eenew(1,l,k,iti)+eenew(2,l,k,iti)*cost1
+ EE(l,k,i-2)=sint1sq*aux
+ gtEE(l,k,i-2)=sint1cost1*aux-sint1cub*eenew(2,l,k,iti)
+ enddo
+ enddo
+ EE(1,1,i-2)=EE(1,1,i-2)+e0new(1,iti)*cost1
+ EE(1,2,i-2)=EE(1,2,i-2)+e0new(2,iti)+e0new(3,iti)*cost1
+ EE(2,1,i-2)=EE(2,1,i-2)+e0new(2,iti)*cost1+e0new(3,iti)
+ EE(2,2,i-2)=EE(2,2,i-2)-e0new(1,iti)
+ gtEE(1,1,i-2)=gtEE(1,1,i-2)-e0new(1,iti)*sint1
+ gtEE(1,2,i-2)=gtEE(1,2,i-2)-e0new(3,iti)*sint1
+ gtEE(2,1,i-2)=gtEE(2,1,i-2)-e0new(2,iti)*sint1
+! print *,"after ee"
+
+!c b1tilde(1,i-2)=b1(1,i-2)
+!c b1tilde(2,i-2)=-b1(2,i-2)
+!c b2tilde(1,i-2)=b2(1,i-2)
+!c b2tilde(2,i-2)=-b2(2,i-2)
+#ifdef DEBUG
+ write (iout,*) 'i=',i-2,gtb1(2,i-2),gtb1(1,i-2)
+ write(iout,*) 'b1=',(b1(k,i-2),k=1,2)
+ write(iout,*) 'b2=',(b2(k,i-2),k=1,2)
+ write (iout,*) 'theta=', theta(i-1)
+#endif
+#else
+ if (i.gt. nnt+2 .and. i.lt.nct+2) then
+! write(iout,*) "i,",molnum(i)
+! print *, "i,",molnum(i),i,itype(i-2,1)
+ if (molnum(i).eq.1) then
+ iti = itype2loc(itype(i-2,1))
+ else
+ iti=nloctyp
+ endif
+ else
+ iti=nloctyp
+ endif
+!c write (iout,*) "i",i-1," itype",itype(i-2)," iti",iti
+!c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
+ if (i.gt. nnt+1 .and. i.lt.nct+1) then
+ iti1 = itype2loc(itype(i-1,1))
+ else
+ iti1=nloctyp
+ endif
+! print *,i,iti
+ b1(1,i-2)=b(3,iti)
+ b1(2,i-2)=b(5,iti)
+ b2(1,i-2)=b(2,iti)
+ b2(2,i-2)=b(4,iti)
+ do k=1,2
+ do l=1,2
+ CC(k,l,i-2)=ccold(k,l,iti)
+ DD(k,l,i-2)=ddold(k,l,iti)
+ EE(k,l,i-2)=eeold(k,l,iti)
+ enddo
+ enddo
+#endif
+ b1tilde(1,i-2)= b1(1,i-2)
+ b1tilde(2,i-2)=-b1(2,i-2)
+ b2tilde(1,i-2)= b2(1,i-2)
+ b2tilde(2,i-2)=-b2(2,i-2)
+!c
+ Ctilde(1,1,i-2)= CC(1,1,i-2)
+ Ctilde(1,2,i-2)= CC(1,2,i-2)
+ Ctilde(2,1,i-2)=-CC(2,1,i-2)
+ Ctilde(2,2,i-2)=-CC(2,2,i-2)
+!c
+ Dtilde(1,1,i-2)= DD(1,1,i-2)
+ Dtilde(1,2,i-2)= DD(1,2,i-2)
+ Dtilde(2,1,i-2)=-DD(2,1,i-2)
+ Dtilde(2,2,i-2)=-DD(2,2,i-2)
+ enddo
#ifdef PARMAT
do i=ivec_start+2,ivec_end+2
#else
do i=3,nres+1
#endif
+
! print *,i,"i"
if (i .lt. nres+1) then
sin1=dsin(phi(i))
if (itype(i-2,1).eq.0) then
iti=ntortyp+1
else
- iti = itortyp(itype(i-2,1))
+ iti = itype2loc(itype(i-2,1))
endif
else
- iti=ntortyp+1
+ iti=nloctyp
endif
! if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
if (i.gt. nnt+1 .and. i.lt.nct+1) then
if (itype(i-1,1).eq.0) then
- iti1=ntortyp+1
+ iti1=nloctyp
else
- iti1 = itortyp(itype(i-1,1))
+ iti1 = itype2loc(itype(i-1,1))
endif
else
- iti1=ntortyp+1
+ iti1=nloctyp
endif
! print *,iti,i,"iti",iti1,itype(i-1,1),itype(i-2,1)
!d write (iout,*) '*******i',i,' iti1',iti
-!d write (iout,*) 'b1',b1(:,iti)
-!d write (iout,*) 'b2',b2(:,iti)
+! write (iout,*) 'b1',b1(:,iti)
+! write (iout,*) 'b2',b2(:,i-2)
!d write (iout,*) 'Ug',Ug(:,:,i-2)
! if (i .gt. iatel_s+2) then
if (i .gt. nnt+2) then
- call matvec2(Ug(1,1,i-2),b2(1,iti),Ub2(1,i-2))
- call matmat2(EE(1,1,iti),Ug(1,1,i-2),EUg(1,1,i-2))
+ call matvec2(Ug(1,1,i-2),b2(1,i-2),Ub2(1,i-2))
+#ifdef NEWCORR
+ call matvec2(Ug(1,1,i-2),gtb2(1,i-2),gUb2(1,i-2))
+!c write (iout,*) Ug(1,1,i-2),gtb2(1,i-2),gUb2(1,i-2),"chuj"
+#endif
+
+ call matmat2(EE(1,1,i-2),Ug(1,1,i-2),EUg(1,1,i-2))
+ call matmat2(gtEE(1,1,i-2),Ug(1,1,i-2),gtEUg(1,1,i-2))
if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0) &
then
- call matmat2(CC(1,1,iti),Ug(1,1,i-2),CUg(1,1,i-2))
- call matmat2(DD(1,1,iti),Ug(1,1,i-2),DUg(1,1,i-2))
- call matmat2(Dtilde(1,1,iti),Ug2(1,1,i-2),DtUg2(1,1,i-2))
- call matvec2(Ctilde(1,1,iti1),obrot(1,i-2),Ctobr(1,i-2))
- call matvec2(Dtilde(1,1,iti),obrot2(1,i-2),Dtobr2(1,i-2))
+ call matmat2(CC(1,1,i-2),Ug(1,1,i-2),CUg(1,1,i-2))
+ call matmat2(DD(1,1,i-2),Ug(1,1,i-2),DUg(1,1,i-2))
+ call matmat2(Dtilde(1,1,i-2),Ug2(1,1,i-2),DtUg2(1,1,i-2))
+ call matvec2(Ctilde(1,1,i-1),obrot(1,i-2),Ctobr(1,i-2))
+ call matvec2(Dtilde(1,1,i-2),obrot2(1,i-2),Dtobr2(1,i-2))
endif
else
do k=1,2
enddo
enddo
endif
- call matvec2(Ugder(1,1,i-2),b2(1,iti),Ub2der(1,i-2))
- call matmat2(EE(1,1,iti),Ugder(1,1,i-2),EUgder(1,1,i-2))
+ call matvec2(Ugder(1,1,i-2),b2(1,i-2),Ub2der(1,i-2))
+ call matmat2(EE(1,1,i-2),Ugder(1,1,i-2),EUgder(1,1,i-2))
do k=1,2
muder(k,i-2)=Ub2der(k,i-2)
enddo
! if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
if (i.gt. nnt+1 .and. i.lt.nct+1) then
if (itype(i-1,1).eq.0) then
- iti1=ntortyp+1
+ iti1=nloctyp
elseif (itype(i-1,1).le.ntyp) then
- iti1 = itortyp(itype(i-1,1))
+ iti1 = itype2loc(itype(i-1,1))
else
- iti1=ntortyp+1
+ iti1=nloctyp
endif
else
- iti1=ntortyp+1
+ iti1=nloctyp
endif
do k=1,2
- mu(k,i-2)=Ub2(k,i-2)+b1(k,iti1)
+ mu(k,i-2)=Ub2(k,i-2)+b1(k,i-1)
enddo
-! if (energy_dec) write (iout,*) 'Ub2 ',i,Ub2(:,i-2)
-! if (energy_dec) write (iout,*) 'b1 ',iti1,b1(:,iti1)
-! if (energy_dec) write (iout,*) 'mu ',i,iti1,mu(:,i-2)
+ if (energy_dec) write (iout,*) 'Ub2 ',i,Ub2(:,i-2)
+ if (energy_dec) write (iout,*) 'b1 ',iti1,b1(:,i-1)
+ if (energy_dec) write (iout,*) 'mu ',i,iti1,mu(:,i-2)
!d write (iout,*) 'mu1',mu1(:,i-2)
!d write (iout,*) 'mu2',mu2(:,i-2)
if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.wcorr6.gt.0.0d0) &
then
- call matmat2(CC(1,1,iti1),Ugder(1,1,i-2),CUgder(1,1,i-2))
- call matmat2(DD(1,1,iti),Ugder(1,1,i-2),DUgder(1,1,i-2))
- call matmat2(Dtilde(1,1,iti),Ug2der(1,1,i-2),DtUg2der(1,1,i-2))
- call matvec2(Ctilde(1,1,iti1),obrot_der(1,i-2),Ctobrder(1,i-2))
- call matvec2(Dtilde(1,1,iti),obrot2_der(1,i-2),Dtobr2der(1,i-2))
+ call matmat2(CC(1,1,i-1),Ugder(1,1,i-2),CUgder(1,1,i-2))
+ call matmat2(DD(1,1,i-2),Ugder(1,1,i-2),DUgder(1,1,i-2))
+ call matmat2(Dtilde(1,1,i-2),Ug2der(1,1,i-2),DtUg2der(1,1,i-2))
+ call matvec2(Ctilde(1,1,i-1),obrot_der(1,i-2),Ctobrder(1,i-2))
+ call matvec2(Dtilde(1,1,i-2),obrot2_der(1,i-2),Dtobr2der(1,i-2))
! Vectors and matrices dependent on a single virtual-bond dihedral.
- call matvec2(DD(1,1,iti),b1tilde(1,iti1),auxvec(1))
+ call matvec2(DD(1,1,i-2),b1tilde(1,i-1),auxvec(1))
call matvec2(Ug2(1,1,i-2),auxvec(1),Ug2Db1t(1,i-2))
call matvec2(Ug2der(1,1,i-2),auxvec(1),Ug2Db1tder(1,i-2))
- call matvec2(CC(1,1,iti1),Ub2(1,i-2),CUgb2(1,i-2))
- call matvec2(CC(1,1,iti1),Ub2der(1,i-2),CUgb2der(1,i-2))
+ call matvec2(CC(1,1,i-1),Ub2(1,i-2),CUgb2(1,i-2))
+ call matvec2(CC(1,1,i-1),Ub2der(1,i-2),CUgb2der(1,i-2))
call matmat2(EUg(1,1,i-2),CC(1,1,iti1),EUgC(1,1,i-2))
call matmat2(EUgder(1,1,i-2),CC(1,1,iti1),EUgCder(1,1,i-2))
call matmat2(EUg(1,1,i-2),DD(1,1,iti1),EUgD(1,1,i-2))
real(kind=8),dimension(2,2) :: acipa !el,a_temp
!el real(kind=8),dimension(3,4) :: agg,aggi,aggi1,aggj,aggj1
real(kind=8),dimension(4) :: muij
+ real(kind=8) :: geel_loc_ij,geel_loc_ji
real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
dist_temp, dist_init,rlocshield,fracinbuf
integer xshift,yshift,zshift,ilist,iresshield
do l=1,2
kkk=kkk+1
muij(kkk)=mu(k,i)*mu(l,j)
+#ifdef NEWCORR
+ gmuij1(kkk)=gtb1(k,i+1)*mu(l,j)
+!c write(iout,*) 'k=',k,i,gtb1(k,i+1),gtb1(k,i+1)*mu(l,j)
+ gmuij2(kkk)=gUb2(k,i)*mu(l,j)
+ gmuji1(kkk)=mu(k,i)*gtb1(l,j+1)
+!c write(iout,*) 'l=',l,j,gtb1(l,j+1),gtb1(l,j+1)*mu(k,i)
+ gmuji2(kkk)=mu(k,i)*gUb2(l,j)
+#endif
+
enddo
enddo
!d write (iout,*) 'EELEC: i',i,' j',j
enddo
endif
+#ifdef NEWCORR
+ geel_loc_ij=(a22*gmuij1(1)&
+ +a23*gmuij1(2)&
+ +a32*gmuij1(3)&
+ +a33*gmuij1(4))&
+ *fac_shield(i)*fac_shield(j)&
+ *sss_ele_cut
+
+!c write(iout,*) "derivative over thatai"
+!c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
+!c & a33*gmuij1(4)
+ gloc(nphi+i,icg)=gloc(nphi+i,icg)+&
+ geel_loc_ij*wel_loc
+!c write(iout,*) "derivative over thatai-1"
+!c write(iout,*) a22*gmuij2(1), a23*gmuij2(2) ,a32*gmuij2(3),
+!c & a33*gmuij2(4)
+ geel_loc_ij=&
+ a22*gmuij2(1)&
+ +a23*gmuij2(2)&
+ +a32*gmuij2(3)&
+ +a33*gmuij2(4)
+ gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+&
+ geel_loc_ij*wel_loc&
+ *fac_shield(i)*fac_shield(j)&
+ *sss_ele_cut
+
+
+!c Derivative over j residue
+ geel_loc_ji=a22*gmuji1(1)&
+ +a23*gmuji1(2)&
+ +a32*gmuji1(3)&
+ +a33*gmuji1(4)
+!c write(iout,*) "derivative over thataj"
+!c write(iout,*) a22*gmuji1(1), a23*gmuji1(2) ,a32*gmuji1(3),
+!c & a33*gmuji1(4)
+
+ gloc(nphi+j,icg)=gloc(nphi+j,icg)+&
+ geel_loc_ji*wel_loc&
+ *fac_shield(i)*fac_shield(j)&
+ *sss_ele_cut
+
+
+ geel_loc_ji=&
+ +a22*gmuji2(1)&
+ +a23*gmuji2(2)&
+ +a32*gmuji2(3)&
+ +a33*gmuji2(4)
+!c write(iout,*) "derivative over thataj-1"
+!c write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3),
+!c & a33*gmuji2(4)
+ gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+&
+ geel_loc_ji*wel_loc&
+ *fac_shield(i)*fac_shield(j)&
+ *sss_ele_cut
+#endif
! write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
! eel_loc_ij=0.0
! include 'COMMON.CONTROL'
real(kind=8),dimension(3) :: ggg
real(kind=8),dimension(2,2) :: auxmat,auxmat1,auxmat2,pizda,&
- e1t,e2t,e3t,e1tder,e2tder,e3tder,e1a,ae3,ae3e2
+ e1t,e2t,e3t,e1tder,e2tder,e3tder,e1a,ae3,ae3e2,gpizda1,&
+ gpizda2,auxgmat1,auxgmatt1,auxgmat2,auxgmatt2
+
real(kind=8),dimension(2) :: auxvec,auxvec1
!el real(kind=8),dimension(3,4) :: agg,aggi,aggi1,aggj,aggj1
real(kind=8),dimension(2,2) :: auxmat3 !el, a_temp
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!d call checkint_turn3(i,a_temp,eello_turn3_num)
call matmat2(EUg(1,1,i+1),EUg(1,1,i+2),auxmat(1,1))
+ call matmat2(gtEUg(1,1,i+1),EUg(1,1,i+2),auxgmat1(1,1))
+ call matmat2(EUg(1,1,i+1),gtEUg(1,1,i+2),auxgmat2(1,1))
call transpose2(auxmat(1,1),auxmat1(1,1))
+ call transpose2(auxgmat1(1,1),auxgmatt1(1,1))
+ call transpose2(auxgmat2(1,1),auxgmatt2(1,1))
call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1))
+ call matmat2(a_temp(1,1),auxgmatt1(1,1),gpizda1(1,1))
+ call matmat2(a_temp(1,1),auxgmatt2(1,1),gpizda2(1,1))
+
if (shield_mode.eq.0) then
fac_shield(i)=1.0d0
fac_shield(j)=1.0d0
if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') &
'eturn3',i,j,0.5d0*(pizda(1,1)+pizda(2,2))
+!C#ifdef NEWCORR
+!C Derivatives in theta
+ gloc(nphi+i,icg)=gloc(nphi+i,icg) &
+ +0.5d0*(gpizda1(1,1)+gpizda1(2,2))*wturn3&
+ *fac_shield(i)*fac_shield(j)
+ gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)&
+ +0.5d0*(gpizda2(1,1)+gpizda2(2,2))*wturn3&
+ *fac_shield(i)*fac_shield(j)
+!C#endif
+
+
+
if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. &
(shield_mode.gt.0)) then
!C print *,i,j
! include 'COMMON.CONTROL'
real(kind=8),dimension(3) :: ggg
real(kind=8),dimension(2,2) :: auxmat,auxmat1,auxmat2,pizda,&
- e1t,e2t,e3t,e1tder,e2tder,e3tder,e1a,ae3,ae3e2
- real(kind=8),dimension(2) :: auxvec,auxvec1
+ e1t,e2t,e3t,e1tder,e2tder,e3tder,e1a,ae3,ae3e2,&
+ gte1t,gte2t,gte3t,&
+ gte1a,gtae3,gtae3e2, ae3gte2,&
+ gtEpizda1,gtEpizda2,gtEpizda3
+
+ real(kind=8),dimension(2) :: auxvec,auxvec1,auxgEvec1,auxgEvec2,&
+ auxgEvec3,auxgvec
+
!el real(kind=8),dimension(3,4) :: agg,aggi,aggi1,aggj,aggj1
real(kind=8),dimension(2,2) :: auxmat3 !el a_temp
!el real(kind=8) :: a22,a23,a32,a33,dxi,dyi,dzi,dx_normi,dy_normi,&
!el local variables
integer :: i,j,iti1,iti2,iti3,l,k,ilist,iresshield
real(kind=8) :: eello_turn4,s1,s2,s3,zj,fracinbuf,eello_t4,&
- rlocshield
+ rlocshield,gs23,gs32,gsE13,gs13,gs21,gsE31,gsEE1,gsEE2,gsEE3
j=i+3
! if (j.ne.20) return
a_temp(1,2)=a23
a_temp(2,1)=a32
a_temp(2,2)=a33
- iti1=itortyp(itype(i+1,1))
- iti2=itortyp(itype(i+2,1))
- iti3=itortyp(itype(i+3,1))
+ iti1=i+1
+ iti2=i+2
+ iti3=i+3
! write(iout,*) "iti1",iti1," iti2",iti2," iti3",iti3
call transpose2(EUg(1,1,i+1),e1t(1,1))
call transpose2(Eug(1,1,i+2),e2t(1,1))
call transpose2(Eug(1,1,i+3),e3t(1,1))
+!C Ematrix derivative in theta
+ call transpose2(gtEUg(1,1,i+1),gte1t(1,1))
+ call transpose2(gtEug(1,1,i+2),gte2t(1,1))
+ call transpose2(gtEug(1,1,i+3),gte3t(1,1))
+
call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1))
call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1))
+ call matmat2(gte1t(1,1),a_temp(1,1),gte1a(1,1))
+ call matvec2(e1a(1,1),gUb2(1,i+3),auxgvec(1))
+!c auxalary matrix of E i+1
+ call matvec2(gte1a(1,1),Ub2(1,i+3),auxgEvec1(1))
s1=scalar2(b1(1,iti2),auxvec(1))
+!c derivative of theta i+2 with constant i+3
+ gs23=scalar2(gtb1(1,i+2),auxvec(1))
+!c derivative of theta i+2 with constant i+2
+ gs32=scalar2(b1(1,i+2),auxgvec(1))
+!c derivative of E matix in theta of i+1
+ gsE13=scalar2(b1(1,i+2),auxgEvec1(1))
+
call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1))
+ call matmat2(a_temp(1,1),gte3t(1,1),gtae3(1,1))
call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1))
- s2=scalar2(b1(1,iti1),auxvec(1))
+!c auxilary matrix auxgvec of Ub2 with constant E matirx
+ call matvec2(ae3(1,1),gUb2(1,i+2),auxgvec(1))
+!c auxilary matrix auxgEvec1 of E matix with Ub2 constant
+ call matvec2(gtae3(1,1),Ub2(1,i+2),auxgEvec3(1))
+ s2=scalar2(b1(1,i+1),auxvec(1))
+!c derivative of theta i+1 with constant i+3
+ gs13=scalar2(gtb1(1,i+1),auxvec(1))
+!c derivative of theta i+2 with constant i+1
+ gs21=scalar2(b1(1,i+1),auxgvec(1))
+!c derivative of theta i+3 with constant i+1
+ gsE31=scalar2(b1(1,i+1),auxgEvec3(1))
+
call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
+ call matmat2(gtae3(1,1),e2t(1,1),gtae3e2(1,1))
+!c ae3gte2 is derivative over i+2
+ call matmat2(ae3(1,1),gte2t(1,1),ae3gte2(1,1))
+
call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
+ call matmat2(ae3e2(1,1),gte1t(1,1),gtEpizda1(1,1))
+!c i+2
+ call matmat2(ae3gte2(1,1),e1t(1,1),gtEpizda2(1,1))
+!c i+3
+ call matmat2(gtae3e2(1,1),e1t(1,1),gtEpizda3(1,1))
+
s3=0.5d0*(pizda(1,1)+pizda(2,2))
+ gsEE1=0.5d0*(gtEpizda1(1,1)+gtEpizda1(2,2))
+ gsEE2=0.5d0*(gtEpizda2(1,1)+gtEpizda2(2,2))
+ gsEE3=0.5d0*(gtEpizda3(1,1)+gtEpizda3(2,2))
if (shield_mode.eq.0) then
fac_shield(i)=1.0
fac_shield(j)=1.0
! print *,"gshieldc_t4(k,j+1)",j,gshieldc_t4(k,j+1)
enddo
endif
-
+#ifdef NEWCORR
+ gloc(nphi+i,icg)=gloc(nphi+i,icg)&
+ -(gs13+gsE13+gsEE1)*wturn4&
+ *fac_shield(i)*fac_shield(j)
+ gloc(nphi+i+1,icg)= gloc(nphi+i+1,icg)&
+ -(gs23+gs21+gsEE2)*wturn4&
+ *fac_shield(i)*fac_shield(j)
+
+ gloc(nphi+i+2,icg)= gloc(nphi+i+2,icg)&
+ -(gs32+gsE31+gsEE3)*wturn4&
+ *fac_shield(i)*fac_shield(j)
+
+!c gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)-
+!c & gs2
+#endif
if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') &
'eturn4',i,j,-(s1+s2+s3)
!d write (2,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3),
call transpose2(EUgder(1,1,i+1),e1tder(1,1))
call matmat2(e1tder(1,1),a_temp(1,1),auxmat(1,1))
call matvec2(auxmat(1,1),Ub2(1,i+3),auxvec(1))
- s1=scalar2(b1(1,iti2),auxvec(1))
+ s1=scalar2(b1(1,i+1),auxvec(1))
call matmat2(ae3e2(1,1),e1tder(1,1),pizda(1,1))
s3=0.5d0*(pizda(1,1)+pizda(2,2))
gel_loc_turn4(i)=gel_loc_turn4(i)-(s1+s3) &
end subroutine theteng
#else
!-----------------------------------------------------------------------------
- subroutine ebend(etheta,ethetacnstr)
+ subroutine ebend(etheta)
!
! Evaluate the virtual-bond-angle energy given the virtual-bond dihedral
! angles gamma and its derivatives in consecutive thetas and gammas.
enddo
!-----------thete constrains
! if (tor_mode.ne.2) then
- ethetacnstr=0.0d0
- print *,ithetaconstr_start,ithetaconstr_end,"TU"
- do i=ithetaconstr_start,ithetaconstr_end
- itheta=itheta_constr(i)
- thetiii=theta(itheta)
- difi=pinorm(thetiii-theta_constr0(i))
- if (difi.gt.theta_drange(i)) then
- difi=difi-theta_drange(i)
- ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4
- gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) &
- +for_thet_constr(i)*difi**3
- else if (difi.lt.-drange(i)) then
- difi=difi+drange(i)
- ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4
- gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) &
- +for_thet_constr(i)*difi**3
- else
- difi=0.0
- endif
- if (energy_dec) then
- write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc", &
- i,itheta,rad2deg*thetiii, &
- rad2deg*theta_constr0(i), rad2deg*theta_drange(i), &
- rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4, &
- gloc(itheta+nphi-2,icg)
- endif
- enddo
-! endif
return
end subroutine ebend
end subroutine etor_d
#else
!-----------------------------------------------------------------------------
- subroutine etor(etors,edihcnstr)
+ subroutine etor(etors)
! implicit real*8 (a-h,o-z)
! include 'DIMENSIONS'
! include 'COMMON.VAR'
! write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg)
enddo
! 6/20/98 - dihedral angle constraints
- edihcnstr=0.0d0
-! do i=1,ndih_constr
+ return
+ end subroutine etor
+!C The rigorous attempt to derive energy function
+!-------------------------------------------------------------------------------------------
+ subroutine etor_kcc(etors)
+ double precision c1(0:maxval_kcc),c2(0:maxval_kcc)
+ real(kind=8) :: etors,glocig,glocit1,glocit2,sinthet1,&
+ sinthet2,costhet1,costhet2,sint1t2,sint1t2n,phii,sinphi,cosphi,&
+ sint1t2n1,sumvalc,gradvalct1,gradvalct2,sumvals,gradvalst1,&
+ gradvalst2,etori
+ logical lprn
+ integer :: i,j,itori,itori1,nval,k,l
+
+ if (lprn) write (iout,*) "etor_kcc tor_mode",tor_mode
+ etors=0.0D0
+ do i=iphi_start,iphi_end
+!C ANY TWO ARE DUMMY ATOMS in row CYCLE
+!c if (((itype(i-3).eq.ntyp1).and.(itype(i-2).eq.ntyp1)).or.
+!c & ((itype(i-2).eq.ntyp1).and.(itype(i-1).eq.ntyp1)) .or.
+!c & ((itype(i-1).eq.ntyp1).and.(itype(i).eq.ntyp1))) cycle
+ if (itype(i-2,1).eq.ntyp1.or. itype(i-1,1).eq.ntyp1 &
+ .or. itype(i,1).eq.ntyp1 .or. itype(i-3,1).eq.ntyp1) cycle
+ itori=itortyp(itype(i-2,1))
+ itori1=itortyp(itype(i-1,1))
+ phii=phi(i)
+ glocig=0.0D0
+ glocit1=0.0d0
+ glocit2=0.0d0
+!C to avoid multiple devision by 2
+!c theti22=0.5d0*theta(i)
+!C theta 12 is the theta_1 /2
+!C theta 22 is theta_2 /2
+!c theti12=0.5d0*theta(i-1)
+!C and appropriate sinus function
+ sinthet1=dsin(theta(i-1))
+ sinthet2=dsin(theta(i))
+ costhet1=dcos(theta(i-1))
+ costhet2=dcos(theta(i))
+!C to speed up lets store its mutliplication
+ sint1t2=sinthet2*sinthet1
+ sint1t2n=1.0d0
+!C \sum_{i=1}^n (sin(theta_1) * sin(theta_2))^n * (c_n* cos(n*gamma)
+!C +d_n*sin(n*gamma)) *
+!C \sum_{i=1}^m (1+a_m*Tb_m(cos(theta_1 /2))+b_m*Tb_m(cos(theta_2 /2)))
+!C we have two sum 1) Non-Chebyshev which is with n and gamma
+ nval=nterm_kcc_Tb(itori,itori1)
+ c1(0)=0.0d0
+ c2(0)=0.0d0
+ c1(1)=1.0d0
+ c2(1)=1.0d0
+ do j=2,nval
+ c1(j)=c1(j-1)*costhet1
+ c2(j)=c2(j-1)*costhet2
+ enddo
+ etori=0.0d0
+
+ do j=1,nterm_kcc(itori,itori1)
+ cosphi=dcos(j*phii)
+ sinphi=dsin(j*phii)
+ sint1t2n1=sint1t2n
+ sint1t2n=sint1t2n*sint1t2
+ sumvalc=0.0d0
+ gradvalct1=0.0d0
+ gradvalct2=0.0d0
+ do k=1,nval
+ do l=1,nval
+ sumvalc=sumvalc+v1_kcc(l,k,j,itori1,itori)*c1(k)*c2(l)
+ gradvalct1=gradvalct1+ &
+ (k-1)*v1_kcc(l,k,j,itori1,itori)*c1(k-1)*c2(l)
+ gradvalct2=gradvalct2+ &
+ (l-1)*v1_kcc(l,k,j,itori1,itori)*c1(k)*c2(l-1)
+ enddo
+ enddo
+ gradvalct1=-gradvalct1*sinthet1
+ gradvalct2=-gradvalct2*sinthet2
+ sumvals=0.0d0
+ gradvalst1=0.0d0
+ gradvalst2=0.0d0
+ do k=1,nval
+ do l=1,nval
+ sumvals=sumvals+v2_kcc(l,k,j,itori1,itori)*c1(k)*c2(l)
+ gradvalst1=gradvalst1+ &
+ (k-1)*v2_kcc(l,k,j,itori1,itori)*c1(k-1)*c2(l)
+ gradvalst2=gradvalst2+ &
+ (l-1)*v2_kcc(l,k,j,itori1,itori)*c1(k)*c2(l-1)
+ enddo
+ enddo
+ gradvalst1=-gradvalst1*sinthet1
+ gradvalst2=-gradvalst2*sinthet2
+ if (lprn) write (iout,*)j,"sumvalc",sumvalc," sumvals",sumvals
+ etori=etori+sint1t2n*(sumvalc*cosphi+sumvals*sinphi)
+!C glocig is the gradient local i site in gamma
+ glocig=glocig+j*sint1t2n*(sumvals*cosphi-sumvalc*sinphi)
+!C now gradient over theta_1
+ glocit1=glocit1+sint1t2n*(gradvalct1*cosphi+gradvalst1*sinphi)&
+ +j*sint1t2n1*costhet1*sinthet2*(sumvalc*cosphi+sumvals*sinphi)
+ glocit2=glocit2+sint1t2n*(gradvalct2*cosphi+gradvalst2*sinphi)&
+ +j*sint1t2n1*sinthet1*costhet2*(sumvalc*cosphi+sumvals*sinphi)
+ enddo ! j
+ etors=etors+etori
+ gloc(i-3,icg)=gloc(i-3,icg)+wtor*glocig
+!C derivative over theta1
+ gloc(nphi+i-3,icg)=gloc(nphi+i-3,icg)+wtor*glocit1
+!C now derivative over theta2
+ gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wtor*glocit2
+ if (lprn) then
+ write (iout,*) i-2,i-1,itype(i-2,1),itype(i-1,1),itori,itori1,&
+ theta(i-1)*rad2deg,theta(i)*rad2deg,phii*rad2deg,etori
+ write (iout,*) "c1",(c1(k),k=0,nval), &
+ " c2",(c2(k),k=0,nval)
+ endif
+ enddo
+ return
+ end subroutine etor_kcc
+!------------------------------------------------------------------------------
+
+ subroutine etor_constr(edihcnstr)
+ real(kind=8) :: etors,edihcnstr
+ logical :: lprn
+!el local variables
+ integer :: i,j,iblock,itori,itori1
+ real(kind=8) :: phii,gloci,v1ij,v2ij,cosphi,sinphi,&
+ vl1ij,vl2ij,vl3ij,pom1,difi,etors_ii,pom,&
+ gaudih_i,gauder_i,s,cos_i,dexpcos_i
+
+ if (raw_psipred) then
+ do i=idihconstr_start,idihconstr_end
+ itori=idih_constr(i)
+ phii=phi(itori)
+ gaudih_i=vpsipred(1,i)
+ gauder_i=0.0d0
+ do j=1,2
+ s = sdihed(j,i)
+ cos_i=(1.0d0-dcos(phii-phibound(j,i)))/s**2
+ dexpcos_i=dexp(-cos_i*cos_i)
+ gaudih_i=gaudih_i+vpsipred(j+1,i)*dexpcos_i
+ gauder_i=gauder_i-2*vpsipred(j+1,i)*dsin(phii-phibound(j,i)) &
+ *cos_i*dexpcos_i/s**2
+ enddo
+ edihcnstr=edihcnstr-wdihc*dlog(gaudih_i)
+ gloc(itori-3,icg)=gloc(itori-3,icg)-wdihc*gauder_i/gaudih_i
+ if (energy_dec) &
+ write (iout,'(2i5,f8.3,f8.5,2(f8.5,2f8.3),f10.5)') &
+ i,itori,phii*rad2deg,vpsipred(1,i),vpsipred(2,i),&
+ phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,vpsipred(3,i),&
+ phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,&
+ -wdihc*dlog(gaudih_i)
+ enddo
+ else
+
do i=idihconstr_start,idihconstr_end
itori=idih_constr(i)
phii=phi(itori)
else
difi=0.0
endif
-!d write (iout,'(2i5,4f8.3,2e14.5)') i,itori,rad2deg*phii,
-!d & rad2deg*phi0(i), rad2deg*drange(i),
-!d & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg)
enddo
-!d write (iout,*) 'edihcnstr',edihcnstr
+
+ endif
+
return
- end subroutine etor
+
+ end subroutine etor_constr
!-----------------------------------------------------------------------------
subroutine etor_d(etors_d)
! 6/23/01 Compute double torsional energy
return
end subroutine etor_d
#endif
+
+ subroutine ebend_kcc(etheta)
+ logical lprn
+ double precision thybt1(maxang_kcc),etheta
+ integer :: i,iti,j,ihelp
+ real (kind=8) :: sinthet,costhet,sumth1thyb,gradthybt1
+!C Set lprn=.true. for debugging
+ lprn=energy_dec
+!c lprn=.true.
+!C print *,"wchodze kcc"
+ if (lprn) write (iout,*) "ebend_kcc tor_mode",tor_mode
+ etheta=0.0D0
+ do i=ithet_start,ithet_end
+!c print *,i,itype(i-1),itype(i),itype(i-2)
+ if ((itype(i-1,1).eq.ntyp1).or.itype(i-2,1).eq.ntyp1 &
+ .or.itype(i,1).eq.ntyp1) cycle
+ iti=iabs(itortyp(itype(i-1,1)))
+ sinthet=dsin(theta(i))
+ costhet=dcos(theta(i))
+ do j=1,nbend_kcc_Tb(iti)
+ thybt1(j)=v1bend_chyb(j,iti)
+ enddo
+ sumth1thyb=v1bend_chyb(0,iti)+ &
+ tschebyshev(1,nbend_kcc_Tb(iti),thybt1(1),costhet)
+ if (lprn) write (iout,*) i-1,itype(i-1,1),iti,theta(i)*rad2deg,&
+ sumth1thyb
+ ihelp=nbend_kcc_Tb(iti)-1
+ gradthybt1=gradtschebyshev(0,ihelp,thybt1(1),costhet)
+ etheta=etheta+sumth1thyb
+!C print *,sumth1thyb,gradthybt1,sinthet*(-0.5d0)
+ gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)-wang*gradthybt1*sinthet
+ enddo
+ return
+ end subroutine ebend_kcc
+!c------------
+!c-------------------------------------------------------------------------------------
+ subroutine etheta_constr(ethetacnstr)
+ real (kind=8) :: ethetacnstr,thetiii,difi
+ integer :: i,itheta
+ ethetacnstr=0.0d0
+!C print *,ithetaconstr_start,ithetaconstr_end,"TU"
+ do i=ithetaconstr_start,ithetaconstr_end
+ itheta=itheta_constr(i)
+ thetiii=theta(itheta)
+ difi=pinorm(thetiii-theta_constr0(i))
+ if (difi.gt.theta_drange(i)) then
+ difi=difi-theta_drange(i)
+ ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4
+ gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) &
+ +for_thet_constr(i)*difi**3
+ else if (difi.lt.-drange(i)) then
+ difi=difi+drange(i)
+ ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4
+ gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) &
+ +for_thet_constr(i)*difi**3
+ else
+ difi=0.0
+ endif
+ if (energy_dec) then
+ write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc",&
+ i,itheta,rad2deg*thetiii,&
+ rad2deg*theta_constr0(i), rad2deg*theta_drange(i),&
+ rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4,&
+ gloc(itheta+nphi-2,icg)
+ endif
+ enddo
+ return
+ end subroutine etheta_constr
+
!-----------------------------------------------------------------------------
subroutine eback_sc_corr(esccor)
! 7/21/2007 Correlations between the backbone-local and side-chain-local
iti1 = itortyp(itype(i+1,1))
if (j.lt.nres-1) then
- itj1 = itortyp(itype(j+1,1))
+ itj1 = itype2loc(itype(j+1,1))
else
- itj1=ntortyp+1
+ itj1=nloctyp
endif
do iii=1,2
dipi(iii,1)=Ub2(iii,i)
#ifdef TIMING
time01=MPI_Wtime()
#endif
+!#define DEBUG
#ifdef DEBUG
write (iout,*) "sum_gradient gvdwc, gvdwx"
do i=1,nres
! write (iout,'(i5,3f10.5,2x,f10.5)')
! & i,(gcorr4_turn(j,i),j=1,3),gel_loc_turn4(i)
! enddo
- write (iout,*) "gvdwc gvdwc_scp gvdwc_scpp"
- do i=1,nres
- write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)') &
- i,(gvdwc(j,i),j=1,3),(gvdwc_scp(j,i),j=1,3),&
- (gvdwc_scpp(j,i),j=1,3)
- enddo
- write (iout,*) "gelc_long gvdwpp gel_loc_long"
- do i=1,nres
- write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)') &
- i,(gelc_long(j,i),j=1,3),(gvdwpp(j,i),j=1,3),&
- (gelc_loc_long(j,i),j=1,3)
- enddo
+! write (iout,*) "gvdwc gvdwc_scp gvdwc_scpp"
+! do i=1,nres
+! write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)') &
+! i,(gvdwc(j,i),j=1,3),(gvdwc_scp(j,i),j=1,3),&
+! (gvdwc_scpp(j,i),j=1,3)
+! enddo
+! write (iout,*) "gelc_long gvdwpp gel_loc_long"
+! do i=1,nres
+! write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)') &
+! i,(gelc_long(j,i),j=1,3),(gvdwpp(j,i),j=1,3),&
+! (gelc_loc_long(j,i),j=1,3)
+! enddo
call flush(iout)
#endif
#ifdef SPLITELE
+wcorr3_nucl*gradcorr3_nucl(j,i) +&
wcatprot* gradpepcat(j,i)+ &
wcatcat*gradcatcat(j,i)+ &
- wscbase*gvdwc_scbase(j,i) &
+ wscbase*gvdwc_scbase(j,i)+ &
wpepbase*gvdwc_pepbase(j,i)+&
wscpho*gvdwc_scpho(j,i)+&
wpeppho*gvdwc_peppho(j,i)
+gradafm(j,i) &
+wliptran*gliptranc(j,i) &
+welec*gshieldc(j,i) &
- +welec*gshieldc_loc(j,) &
+ +welec*gshieldc_loc(j,i) &
+wcorr*gshieldc_ec(j,i) &
+wcorr*gshieldc_loc_ec(j,i) &
+wturn3*gshieldc_t3(j,i) &
enddo
return
end subroutine sc_grad
-#ifdef CRYST_THETA
-!-----------------------------------------------------------------------------
- subroutine mixder(thetai,thet_pred_mean,theta0i,E_tc_t)
- use comm_calcthet
+ subroutine sc_grad_cat
! implicit real*8 (a-h,o-z)
+ use calc_data
! include 'DIMENSIONS'
-! include 'COMMON.LOCAL'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.CALC'
! include 'COMMON.IOUNITS'
-!el real(kind=8) :: term1,term2,termm,diffak,ratak,&
-!el ak,aktc,termpre,termexp,sigc,sig0i,time11,time12,sigcsq,&
-!el delthe0,sig0inv,sigtc,sigsqtc,delthec,
- real(kind=8) :: thetai,thet_pred_mean,theta0i,E_tc_t
- real(kind=8) :: t3,t6,t9,t12,t14,t16,t21,t23,t26,t27,t32,t40
-!el integer :: it
-!el common /calcthet/ term1,term2,termm,diffak,ratak,&
-!el ak,aktc,termpre,termexp,sigc,sig0i,time11,time12,sigcsq,&
-!el delthe0,sig0inv,sigtc,sigsqtc,delthec,it
-!el local variables
+ real(kind=8), dimension(3) :: dcosom1,dcosom2
+! print *,"wchodze"
+ eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1 &
+ +dCAVdOM1+ dGCLdOM1+ dPOLdOM1
+ eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2 &
+ +dCAVdOM2+ dGCLdOM2+ dPOLdOM2
- delthec=thetai-thet_pred_mean
- delthe0=thetai-theta0i
-! "Thank you" to MAPLE (probably spared one day of hand-differentiation).
- t3 = thetai-thet_pred_mean
- t6 = t3**2
- t9 = term1
- t12 = t3*sigcsq
- t14 = t12+t6*sigsqtc
- t16 = 1.0d0
- t21 = thetai-theta0i
- t23 = t21**2
- t26 = term2
- t27 = t21*t26
- t32 = termexp
- t40 = t32**2
- E_tc_t = -((sigcsq+2.D0*t3*sigsqtc)*t9-t14*sigcsq*t3*t16*t9 &
- -aktc*sig0inv*t27)/t32+(t14*t9+aktc*t26)/t40 &
- *(-t12*t9-ak*sig0inv*t27)
- return
- end subroutine mixder
-#endif
-!-----------------------------------------------------------------------------
-! cartder.F
-!-----------------------------------------------------------------------------
- subroutine cartder
+ eom12=evdwij*eps1_om12+eps2der*eps2rt_om12 &
+ -2.0D0*alf12*eps3der+sigder*sigsq_om12&
+ +dCAVdOM12+ dGCLdOM12
+! diagnostics only
+! eom1=0.0d0
+! eom2=0.0d0
+! eom12=evdwij*eps1_om12
+! end diagnostics
+! write (iout,*) "eps2der",eps2der," eps3der",eps3der,&
+! " sigder",sigder
+! write (iout,*) "eps1_om12",eps1_om12," eps2rt_om12",eps2rt_om12
+! write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12
+!C print *,sss_ele_cut,'in sc_grad'
+
+ do k=1,3
+ dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k))
+ dcosom2(k)=rij*(dc_norm(k,j)-om2*erij(k))
+ enddo
+ do k=1,3
+ gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))
+!C print *,'gg',k,gg(k)
+ enddo
+! print *,i,j,gg_lipi(3),gg_lipj(3),sss_ele_cut
+! write (iout,*) "gg",(gg(k),k=1,3)
+ do k=1,3
+ gvdwx(k,i)=gvdwx(k,i)-gg(k) +gg_lipi(k)&
+ +(eom12*(dc_norm(k,j)-om12*dc_norm(k,nres+i)) &
+ +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
+
+ gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)&
+ +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,j)) &
+ +eom2*(erij(k)-om2*dc_norm(k,j)))*dscj_inv
+
+! write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) &
+! +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
+! write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) &
+! +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
+ enddo
+!
+! Calculate the components of the gradient in DC and X
+!
+!grad do k=i,j-1
+!grad do l=1,3
+!grad gvdwc(l,k)=gvdwc(l,k)+gg(l)
+!grad enddo
+!grad enddo
+ do l=1,3
+ gvdwc(l,i)=gvdwc(l,i)-gg(l)
+ gvdwc(l,j)=gvdwc(l,j)+gg(l)
+ enddo
+ end subroutine sc_grad_cat
+
+
+#ifdef CRYST_THETA
+!-----------------------------------------------------------------------------
+ subroutine mixder(thetai,thet_pred_mean,theta0i,E_tc_t)
+
+ use comm_calcthet
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.LOCAL'
+! include 'COMMON.IOUNITS'
+!el real(kind=8) :: term1,term2,termm,diffak,ratak,&
+!el ak,aktc,termpre,termexp,sigc,sig0i,time11,time12,sigcsq,&
+!el delthe0,sig0inv,sigtc,sigsqtc,delthec,
+ real(kind=8) :: thetai,thet_pred_mean,theta0i,E_tc_t
+ real(kind=8) :: t3,t6,t9,t12,t14,t16,t21,t23,t26,t27,t32,t40
+!el integer :: it
+!el common /calcthet/ term1,term2,termm,diffak,ratak,&
+!el ak,aktc,termpre,termexp,sigc,sig0i,time11,time12,sigcsq,&
+!el delthe0,sig0inv,sigtc,sigsqtc,delthec,it
+!el local variables
+
+ delthec=thetai-thet_pred_mean
+ delthe0=thetai-theta0i
+! "Thank you" to MAPLE (probably spared one day of hand-differentiation).
+ t3 = thetai-thet_pred_mean
+ t6 = t3**2
+ t9 = term1
+ t12 = t3*sigcsq
+ t14 = t12+t6*sigsqtc
+ t16 = 1.0d0
+ t21 = thetai-theta0i
+ t23 = t21**2
+ t26 = term2
+ t27 = t21*t26
+ t32 = termexp
+ t40 = t32**2
+ E_tc_t = -((sigcsq+2.D0*t3*sigsqtc)*t9-t14*sigcsq*t3*t16*t9 &
+ -aktc*sig0inv*t27)/t32+(t14*t9+aktc*t26)/t40 &
+ *(-t12*t9-ak*sig0inv*t27)
+ return
+ end subroutine mixder
+#endif
+!-----------------------------------------------------------------------------
+! cartder.F
+!-----------------------------------------------------------------------------
+ subroutine cartder
!-----------------------------------------------------------------------------
! This subroutine calculates the derivatives of the consecutive virtual
! bond vectors and the SC vectors in the virtual-bond angles theta and
! call intcartderiv
! call checkintcartgrad
call zerograd
- aincr=2.0D-5
+ aincr=1.0D-7
write(iout,*) 'Calling CHECK_ECARTINT.',aincr
nf=0
icall=0
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
rij=dsqrt(rrij)
sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
- sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj)))
- sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj)))
+ sss_ele_cut=sscale_ele(1.0d0/(rij))
+ sss_ele_grad=sscagrad_ele(1.0d0/(rij))
sss_grad=sscale_grad(1.0d0/(rij*sigmaii(itypi,itypj)))
if (sss_ele_cut.le.0.0) cycle
if (sss.lt.1.0d0) then
sigder=fac*sigder
fac=rij*fac
fac=fac+evdwij*(sss_ele_grad/sss_ele_cut&
- /sigma(itypi,itypj)*rij-sss_grad/(1.0-sss)*rij &
+ *rij-sss_grad/(1.0-sss)*rij &
/sigmaii(itypi,itypj))
! fac=0.0d0
! Calculate the radial part of the gradient
rij=dsqrt(rrij)
sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
sss_grad=sscale_grad(1.0d0/(rij*sigmaii(itypi,itypj)))
- sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj)))
- sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj)))
+ sss_ele_cut=sscale_ele(1.0d0/(rij))
+ sss_ele_grad=sscagrad_ele(1.0d0/(rij))
if (sss_ele_cut.le.0.0) cycle
if (sss.gt.0.0d0) then
sigder=fac*sigder
fac=rij*fac
fac=fac+evdwij*(sss_ele_grad/sss_ele_cut&
- /sigma(itypi,itypj)*rij+sss_grad/sss*rij &
+ *rij+sss_grad/sss*rij &
/sigmaii(itypi,itypj))
! fac=0.0d0
!
! Calculate the virtual-bond-angle energy.
!
- call ebend(ebe,ethetacnstr)
-!
! Calculate the SC local energy.
!
call vec_and_deriv
call esc(escloc)
!
+ if (wang.gt.0d0) then
+ if (tor_mode.eq.0) then
+ call ebend(ebe)
+ else
+!C ebend kcc is Kubo cumulant clustered rigorous attemp to derive the
+!C energy function
+ call ebend_kcc(ebe)
+ endif
+ else
+ ebe=0.0d0
+ endif
+ ethetacnstr=0.0d0
+ if (with_theta_constr) call etheta_constr(ethetacnstr)
+
+! write(iout,*) "in etotal afer ebe",ipot
+
+! print *,"Processor",myrank," computed UB"
+!
+! Calculate the SC local energy.
+!
+ call esc(escloc)
+!elwrite(iout,*) "in etotal afer esc",ipot
+! print *,"Processor",myrank," computed USC"
+!
+! Calculate the virtual-bond torsional energy.
+!
+!d print *,'nterm=',nterm
+! if (wtor.gt.0) then
+! call etor(etors,edihcnstr)
+! else
+! etors=0
+! edihcnstr=0
+! endif
+ if (wtor.gt.0.0d0) then
+ if (tor_mode.eq.0) then
+ call etor(etors)
+ else
+!C etor kcc is Kubo cumulant clustered rigorous attemp to derive the
+!C energy function
+ call etor_kcc(etors)
+ endif
+ else
+ etors=0.0d0
+ endif
+ edihcnstr=0.0d0
+ if (ndih_constr.gt.0) call etor_constr(edihcnstr)
+
! Calculate the virtual-bond torsional energy.
!
- call etor(etors,edihcnstr)
!
! 6/23/01 Calculate double-torsional energy
!
+ if ((wtor_d.gt.0.0d0).and.(tor_mode.eq.0)) then
call etor_d(etors_d)
+ endif
!
! 21/5/07 Calculate local sicdechain correlation energy
!
integer :: i,j
if(nres.lt.100) then
- maxconts=nres
+ maxconts=10*nres
elseif(nres.lt.200) then
- maxconts=0.8*nres ! Max. number of contacts per residue
+ maxconts=10*nres ! Max. number of contacts per residue
else
- maxconts=0.6*nres ! (maxconts=maxres/4)
+ maxconts=10*nres ! (maxconts=maxres/4)
endif
maxcont=12*nres ! Max. number of SC contacts
maxvar=6*nres ! Max. number of variables
allocate(Ug2DtEUgder(2,2,2,nres))
allocate(DtUg2EUgder(2,2,2,nres))
!(2,2,2,maxres)
+ allocate(b1(2,nres)) !(2,-maxtor:maxtor)
+ allocate(b2(2,nres)) !(2,-maxtor:maxtor)
+ allocate(b1tilde(2,nres)) !(2,-maxtor:maxtor)
+ allocate(b2tilde(2,nres)) !(2,-maxtor:maxtor)
+
+ allocate(ctilde(2,2,nres))
+ allocate(dtilde(2,2,nres)) !(2,2,-maxtor:maxtor)
+ allocate(gtb1(2,nres))
+ allocate(gtb2(2,nres))
+ allocate(cc(2,2,nres))
+ allocate(dd(2,2,nres))
+ allocate(ee(2,2,nres))
+ allocate(gtcc(2,2,nres))
+ allocate(gtdd(2,2,nres))
+ allocate(gtee(2,2,nres))
+ allocate(gUb2(2,nres))
+ allocate(gteUg(2,2,nres))
+
! common /rotat_old/
allocate(costab(nres))
allocate(sintab(nres))
enddo
! IF ( (wcorr_nucl.gt.0.0d0.or.wcorr3_nucl.gt.0.0d0) .and.
IF ( j.gt.i+1 .and.&
- num_conti.le.maxconts) THEN
+ num_conti.le.maxcont) THEN
!C
!C Calculate the contact function. The ith column of the array JCONT will
!C contain the numbers of atoms that make contacts with the atom I (of numbers
!C greater than I). The arrays FACONT and GACONT will contain the values of
!C the contact function and its derivative.
- r0ij=2.20D0*sigma(itypi,itypj)
+ r0ij=2.20D0*sigma_nucl(itypi,itypj)
!c write (2,*) "ij",i,j," rij",1.0d0/rij," r0ij",r0ij
call gcont(rij,r0ij,1.0D0,0.2d0/r0ij,fcont,fprimcont)
!c write (2,*) "fcont",fcont
if (num_conti.gt.maxconts) then
write (iout,*) 'WARNING - max. # of contacts exceeded;',&
- ' will skip next contacts for this conf.'
+ ' will skip next contacts for this conf.',maxconts
else
jcont_hb(num_conti,i)=j
!c write (iout,*) "num_conti",num_conti,
return
end subroutine ecatcat
!---------------------------------------------------------------------------
- subroutine ecat_prot(ecation_prot)
- integer i,j,k,subchap,itmp,inum
- real(kind=8) :: xi,yi,zi,xj,yj,zj,ract,rcat0,epscalc,r06,r012,&
- r7,r4,ecationcation
- real(kind=8) xj_temp,yj_temp,zj_temp,xj_safe,yj_safe,zj_safe, &
- dist_init,dist_temp,ecation_prot,rcal,rocal, &
- Evan1,Evan2,EC,cm1mag,DASGL,delta,r0p,Epepcat, &
- catl,cml,calpl, Etotal_p, Etotal_m,rtab,wdip,wmodquad,wquad1, &
- wquad2,wvan1,E1,E2,wconst,wvan2,rcpm,dcmag,sin2thet,sinthet, &
- costhet,v1m,v2m,wh2o,wc,rsecp,Ir,Irsecp,Irthrp,Irfourp,Irfiftp,&
- Irsistp,Irseven,Irtwelv,Irthir,dE1dr,dE2dr,dEdcos,wquad2p,opt, &
- rs,rthrp,rfourp,rsixp,reight,Irsixp,Ireight,Irtw,Irfourt, &
- opt1,opt2,opt3,opt4,opt5,opt6,opt7,opt8,opt9,opt10,opt11,opt12,&
- opt13,opt14,opt15,opt16,opt17,opt18,opt19, &
- Equad1,Equad2,dscmag,v1dpv2,dscmag3,constA,constB,Edip
- real(kind=8),dimension(3) ::dEvan1Cmcat,dEvan2Cmcat,dEeleccat,&
- gg,r,EtotalCat,dEtotalCm,dEtotalCalp,dEvan1Cm,dEvan2Cm, &
- dEtotalpep,dEtotalcat_num,dEddci,dEtotalcm_num,dEtotalcalp_num, &
- tab1,tab2,tab3,diff,cm1,sc,p,tcat,talp,cm,drcp,drcp_norm,vcat, &
- v1,v2,v3,myd_norm,dx,vcm,valpha,drdpep,dcosdpep,dcosddci,dEdpep,&
- dEcCat,dEdipCm,dEdipCalp,dEquad1Cat,dEquad1Cm,dEquad1Calp, &
- dEquad2Cat,dEquad2Cm,dEquad2Calpd,Evan1Cat,dEvan1Calp,dEvan2Cat,&
- dEvan2Calp,dEtotalCat,dscvec,dEcCm,dEcCalp,dEdipCat,dEquad2Calp,&
- dEvan1Cat
- real(kind=8),dimension(6) :: vcatprm
- ecation_prot=0.0d0
-! first lets calculate interaction with peptide groups
- if (nres_molec(5).eq.0) return
- wconst=78
- wdip =1.092777950857032D2
- wdip=wdip/wconst
- wmodquad=-2.174122713004870D4
- wmodquad=wmodquad/wconst
- wquad1 = 3.901232068562804D1
- wquad1=wquad1/wconst
- wquad2 = 3
- wquad2=wquad2/wconst
- wvan1 = 0.1
- wvan2 = 6
+! new for K+
+ subroutine ecats_prot_amber(ecations_prot_amber)
+! subroutine ecat_prot2(ecation_prot)
+ use calc_data
+ use comm_momo
+
+ logical :: lprn
+!el local variables
+ integer :: iint,itypi1,subchap,isel,itmp
+ real(kind=8) :: rrij,xi,yi,zi,sig,rij_shift,e1,e2,sigm,epsi
+ real(kind=8) :: evdw
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init,ssgradlipi,ssgradlipj, &
+ sslipi,sslipj,faclip,alpha_sco
+ integer :: ii
+ real(kind=8) :: fracinbuf
+ real (kind=8) :: escpho
+ real (kind=8),dimension(4):: ener
+ real(kind=8) :: b1,b2,egb
+ real(kind=8) :: Fisocav,ECL,Elj,Equad,Epol,eheadtail,&
+ Lambf,&
+ Chif,ChiLambf,Fcav,dFdR,dFdOM1,&
+ ecations_prot_amber,dFdOM2,dFdL,dFdOM12,&
+ federmaus,&
+ d1i,d1j
+! real(kind=8),dimension(3,2)::erhead_tail
+! real(kind=8),dimension(3) :: Rhead_distance,ertail,erhead,Rtail_distance
+ real(kind=8) :: facd4, adler, Fgb, facd3
+ integer troll,jj,istate
+ real (kind=8) :: dcosom1(3),dcosom2(3)
+
+ ecations_prot_amber=0.0D0
+ if (nres_molec(5).eq.0) return
+ eps_out=80.0d0
+! sss_ele_cut=1.0d0
+
itmp=0
do i=1,4
itmp=itmp+nres_molec(i)
enddo
! do i=1,nres_molec(1)-1 ! loop over all peptide groups needs parralelization
do i=ibond_start,ibond_end
-! cycle
- if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1)) cycle ! leave dummy atoms
- xi=0.5d0*(c(1,i)+c(1,i+1))
- yi=0.5d0*(c(2,i)+c(2,i+1))
- zi=0.5d0*(c(3,i)+c(3,i+1))
- xi=mod(xi,boxxsize)
+
+! print *,"I am in EVDW",i
+ itypi=iabs(itype(i,1))
+! if (i.ne.47) cycle
+ if (itypi.eq.ntyp1) cycle
+ itypi1=iabs(itype(i+1,1))
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ xi=dmod(xi,boxxsize)
if (xi.lt.0) xi=xi+boxxsize
- yi=mod(yi,boxysize)
+ yi=dmod(yi,boxysize)
if (yi.lt.0) yi=yi+boxysize
- zi=mod(zi,boxzsize)
+ zi=dmod(zi,boxzsize)
if (zi.lt.0) zi=zi+boxzsize
-
+ dxi=dc_norm(1,nres+i)
+ dyi=dc_norm(2,nres+i)
+ dzi=dc_norm(3,nres+i)
+ dsci_inv=vbld_inv(i+nres)
do j=itmp+1,itmp+nres_molec(5)
+
+! Calculate SC interaction energy.
+ itypj=iabs(itype(j,1))
+ if ((itypj.eq.ntyp1)) cycle
+ CALL elgrad_init_cat(eheadtail,Egb,Ecl,Elj,Equad,Epol)
+
+ dscj_inv=vbld_inv(j+nres)
xj=c(1,j)
yj=c(2,j)
zj=c(3,j)
- xj=dmod(xj,boxxsize)
- if (xj.lt.0) xj=xj+boxxsize
- yj=dmod(yj,boxysize)
- if (yj.lt.0) yj=yj+boxysize
- zj=dmod(zj,boxzsize)
- if (zj.lt.0) zj=zj+boxzsize
- dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- xj_safe=xj
- yj_safe=yj
- zj_safe=zj
- subchap=0
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
+ xj=dmod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=dmod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=dmod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
xj=xj_safe+xshift*boxxsize
yj=yj_safe+yshift*boxysize
zj=zj_safe+zshift*boxzsize
zj_temp=zj
subchap=1
endif
- enddo
- enddo
- enddo
- if (subchap.eq.1) then
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
xj=xj_temp-xi
yj=yj_temp-yi
zj=zj_temp-zi
- else
+ else
xj=xj_safe-xi
yj=yj_safe-yi
zj=zj_safe-zi
- endif
-! enddo
-! enddo
- rcpm = sqrt(xj**2+yj**2+zj**2)
- drcp_norm(1)=xj/rcpm
- drcp_norm(2)=yj/rcpm
- drcp_norm(3)=zj/rcpm
- dcmag=0.0
- do k=1,3
- dcmag=dcmag+dc(k,i)**2
- enddo
- dcmag=dsqrt(dcmag)
+ endif
+
+! dxj = dc_norm( 1, nres+j )
+! dyj = dc_norm( 2, nres+j )
+! dzj = dc_norm( 3, nres+j )
+
+ itypi = itype(i,1)
+ itypj = itype(j,5)
+! Parameters from fitting the analitical expressions to the PMF obtained by umbrella
+! sampling performed with amber package
+! alf1 = 0.0d0
+! alf2 = 0.0d0
+! alf12 = 0.0d0
+! a12sq = rborn(itypi,itypj) * rborn(itypj,itypi)
+ chi1 = chicat(itypi,itypj)
+ chis1 = chiscat(itypi,itypj)
+ chip1 = chippcat(itypi,itypj)
+! chis2 = chis(itypj,itypi)
+! chis12 = chis1 * chis2
+ sig1 = sigmap1cat(itypi,itypj)
+! sig2 = sigmap2(itypi,itypj)
+! alpha factors from Fcav/Gcav
+ b1cav = alphasurcat(1,itypi,itypj)
+ b2cav = alphasurcat(2,itypi,itypj)
+ b3cav = alphasurcat(3,itypi,itypj)
+ b4cav = alphasurcat(4,itypi,itypj)
+
+! used to determine whether we want to do quadrupole calculations
+ eps_in = epsintabcat(itypi,itypj)
+ if (eps_in.eq.0.0) eps_in=1.0
+
+ eps_inout_fac = ( (1.0d0/eps_in) - (1.0d0/eps_out))
+! Rtail = 0.0d0
+
+ DO k = 1, 3
+ ctail(k,1)=c(k,i+nres)
+ ctail(k,2)=c(k,j)
+ END DO
+!c! tail distances will be themselves usefull elswhere
+!c1 (in Gcav, for example)
+ Rtail_distance(1) = ctail( 1, 2 ) - ctail( 1,1 )
+ Rtail_distance(2) = ctail( 2, 2 ) - ctail( 2,1 )
+ Rtail_distance(3) = ctail( 3, 2 ) - ctail( 3,1 )
+ Rtail = dsqrt( &
+ (Rtail_distance(1)*Rtail_distance(1)) &
+ + (Rtail_distance(2)*Rtail_distance(2)) &
+ + (Rtail_distance(3)*Rtail_distance(3)))
+! tail location and distance calculations
+! dhead1
+ d1 = dheadcat(1, 1, itypi, itypj)
+! d2 = dhead(2, 1, itypi, itypj)
+ DO k = 1,3
+! location of polar head is computed by taking hydrophobic centre
+! and moving by a d1 * dc_norm vector
+! see unres publications for very informative images
+ chead(k,1) = c(k, i+nres) + d1 * dc_norm(k, i+nres)
+ chead(k,2) = c(k, j)
+! distance
+! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres))
+! Rsc(k) = Rsc_distance(k) * Rsc_distance(k)
+ Rhead_distance(k) = chead(k,2) - chead(k,1)
+ END DO
+! pitagoras (root of sum of squares)
+ Rhead = dsqrt( &
+ (Rhead_distance(1)*Rhead_distance(1)) &
+ + (Rhead_distance(2)*Rhead_distance(2)) &
+ + (Rhead_distance(3)*Rhead_distance(3)))
+!-------------------------------------------------------------------
+! zero everything that should be zero'ed
+ evdwij = 0.0d0
+ ECL = 0.0d0
+ Elj = 0.0d0
+ Equad = 0.0d0
+ Epol = 0.0d0
+ Fcav=0.0d0
+ eheadtail = 0.0d0
+ dGCLdOM1 = 0.0d0
+ dGCLdOM2 = 0.0d0
+ dGCLdOM12 = 0.0d0
+ dPOLdOM1 = 0.0d0
+ dPOLdOM2 = 0.0d0
+ Fcav = 0.0d0
+ dFdR = 0.0d0
+ dCAVdOM1 = 0.0d0
+ dCAVdOM2 = 0.0d0
+ dCAVdOM12 = 0.0d0
+ dscj_inv = vbld_inv(j+nres)
+! print *,i,j,dscj_inv,dsci_inv
+! rij holds 1/(distance of Calpha atoms)
+ rrij = 1.0D0 / ( xj*xj + yj*yj + zj*zj)
+ rij = dsqrt(rrij)
+ CALL sc_angular
+! this should be in elgrad_init but om's are calculated by sc_angular
+! which in turn is used by older potentials
+! om = omega, sqom = om^2
+ sqom1 = om1 * om1
+ sqom2 = om2 * om2
+ sqom12 = om12 * om12
+
+! now we calculate EGB - Gey-Berne
+! It will be summed up in evdwij and saved in evdw
+ sigsq = 1.0D0 / sigsq
+ sig = sig0ij * dsqrt(sigsq)
+! rij_shift = 1.0D0 / rij - sig + sig0ij
+ rij_shift = Rtail - sig + sig0ij
+ IF (rij_shift.le.0.0D0) THEN
+ evdw = 1.0D20
+ RETURN
+ END IF
+ sigder = -sig * sigsq
+ rij_shift = 1.0D0 / rij_shift
+ fac = rij_shift**expon
+ c1 = fac * fac * aa_aq(itypi,itypj)
+! print *,"ADAM",aa_aq(itypi,itypj)
+
+! c1 = 0.0d0
+ c2 = fac * bb_aq(itypi,itypj)
+! c2 = 0.0d0
+ evdwij = eps1 * eps2rt * eps3rt * ( c1 + c2 )
+ eps2der = eps3rt * evdwij
+ eps3der = eps2rt * evdwij
+! evdwij = 4.0d0 * eps2rt * eps3rt * evdwij
+ evdwij = eps2rt * eps3rt * evdwij
+!#ifdef TSCSC
+! IF (bb_aq(itypi,itypj).gt.0) THEN
+! evdw_p = evdw_p + evdwij
+! ELSE
+! evdw_m = evdw_m + evdwij
+! END IF
+!#else
+ evdw = evdw &
+ + evdwij
+!#endif
+ c1 = c1 * eps1 * eps2rt**2 * eps3rt**2
+ fac = -expon * (c1 + evdwij) * rij_shift
+ sigder = fac * sigder
+! Calculate distance derivative
+ gg(1) = fac
+ gg(2) = fac
+ gg(3) = fac
+
+ fac = chis1 * sqom1 + chis2 * sqom2 &
+ - 2.0d0 * chis12 * om1 * om2 * om12
+ pom = 1.0d0 - chis1 * chis2 * sqom12
+ Lambf = (1.0d0 - (fac / pom))
+ Lambf = dsqrt(Lambf)
+ sparrow = 1.0d0 / dsqrt(sig1**2.0d0 + sig2**2.0d0)
+ Chif = Rtail * sparrow
+ ChiLambf = Chif * Lambf
+ eagle = dsqrt(ChiLambf)
+ bat = ChiLambf ** 11.0d0
+ top = b1cav * ( eagle + b2cav * ChiLambf - b3cav )
+ bot = 1.0d0 + b4cav * (ChiLambf ** 12.0d0)
+ botsq = bot * bot
+ Fcav = top / bot
+
+ dtop = b1cav * ((Lambf / (2.0d0 * eagle)) + (b2cav * Lambf))
+ dbot = 12.0d0 * b4cav * bat * Lambf
+ dFdR = ((dtop * bot - top * dbot) / botsq) * sparrow
+
+ dtop = b1cav * ((Chif / (2.0d0 * eagle)) + (b2cav * Chif))
+ dbot = 12.0d0 * b4cav * bat * Chif
+ eagle = Lambf * pom
+ dFdOM1 = -(chis1 * om1 - chis12 * om2 * om12) / (eagle)
+ dFdOM2 = -(chis2 * om2 - chis12 * om1 * om12) / (eagle)
+ dFdOM12 = chis12 * (chis1 * om1 * om12 - om2) &
+ * (chis2 * om2 * om12 - om1) / (eagle * pom)
+
+ dFdL = ((dtop * bot - top * dbot) / botsq)
+ dCAVdOM1 = dFdL * ( dFdOM1 )
+ dCAVdOM2 = dFdL * ( dFdOM2 )
+ dCAVdOM12 = dFdL * ( dFdOM12 )
+
+ DO k= 1, 3
+ ertail(k) = Rtail_distance(k)/Rtail
+ END DO
+ erdxi = scalar( ertail(1), dC_norm(1,i+nres) )
+ erdxj = scalar( ertail(1), dC_norm(1,j) )
+ facd1 = dtail(1,itypi,itypj) * vbld_inv(i+nres)
+ facd2 = dtail(2,itypi,itypj) * vbld_inv(j+nres)
+ DO k = 1, 3
+ pom = ertail(k)-facd1*(ertail(k)-erdxi*dC_norm(k,i+nres))
+ gvdwx(k,i) = gvdwx(k,i) &
+ - (( dFdR + gg(k) ) * pom)
+ pom = ertail(k)-facd2*(ertail(k)-erdxj*dC_norm(k,j+nres))
+! gvdwx(k,j) = gvdwx(k,j) &
+! + (( dFdR + gg(k) ) * pom)
+ gvdwc(k,i) = gvdwc(k,i) &
+ - (( dFdR + gg(k) ) * ertail(k))
+ gvdwc(k,j) = gvdwc(k,j) &
+ + (( dFdR + gg(k) ) * ertail(k))
+ gg(k) = 0.0d0
+
+!c! Compute head-head and head-tail energies for each state
+ isel = iabs(Qi) + iabs(Qj)
+ IF (isel.eq.0) THEN
+!c! No charges - do nothing
+ eheadtail = 0.0d0
+
+ ELSE IF (isel.eq.1 .and. iabs(Qj).eq.1) THEN
+!c! Nonpolar-charge interactions
+ if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
+ Qi=Qi*2
+ Qij=Qij*2
+ endif
+ if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then
+ Qj=Qj*2
+ Qij=Qij*2
+ endif
+
+ CALL enq_cat(epol)
+ eheadtail = epol
+! eheadtail = 0.0d0
+
+ ELSE IF (isel.eq.3 .and. icharge(itypj).eq.2) THEN
+!c! Dipole-charge interactions
+ if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
+ Qi=Qi*2
+ Qij=Qij*2
+ endif
+ if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then
+ Qj=Qj*2
+ Qij=Qij*2
+ endif
+ CALL edq_cat(ecl, elj, epol)
+ eheadtail = ECL + elj + epol
+! eheadtail = 0.0d0
+
+ ELSE IF ((isel.eq.2.and. &
+ iabs(Qi).eq.1).and. &
+ nstate(itypi,itypj).eq.1) THEN
+
+!c! Same charge-charge interaction ( +/+ or -/- )
+ if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
+ Qi=Qi*2
+ Qij=Qij*2
+ endif
+ if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then
+ Qj=Qj*2
+ Qij=Qij*2
+ endif
+
+ CALL eqq_cat(Ecl,Egb,Epol,Fisocav,Elj)
+ eheadtail = ECL + Egb + Epol + Fisocav + Elj
+! eheadtail = 0.0d0
+
+! ELSE IF ((isel.eq.2.and. &
+! iabs(Qi).eq.1).and. &
+! nstate(itypi,itypj).ne.1) THEN
+!c! Different charge-charge interaction ( +/- or -/+ )
+! if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
+! Qi=Qi*2
+! Qij=Qij*2
+! endif
+! if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then
+! Qj=Qj*2
+! Qij=Qij*2
+! endif
+!
+! CALL energy_quad(istate,eheadtail,Ecl,Egb,Epol,Fisocav,Elj,Equad)
+ END IF ! this endif ends the "catch the gly-gly" at the beggining of Fcav
+ evdw = evdw + Fcav + eheadtail
+
+ IF (energy_dec) write (iout,'(2(1x,a3,i3),3f6.2,10f16.7)') &
+ restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,&
+ 1.0d0/rij,Rtail,Rhead,evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,&
+ Equad,evdwij+Fcav+eheadtail,evdw
+! evdw = evdw + Fcav + eheadtail
+
+! iF (nstate(itypi,itypj).eq.1) THEN
+ CALL sc_grad_cat
+! END IF
+!c!-------------------------------------------------------------------
+!c! NAPISY KONCOWE
+ END DO ! j
+ END DO ! iint
+ END DO ! i
+!c write (iout,*) "Number of loop steps in EGB:",ind
+!c energy_dec=.false.
+! print *,"EVDW KURW",evdw,nres
+
+ return
+ end subroutine ecats_prot_amber
+
+!---------------------------------------------------------------------------
+! old for Ca2+
+ subroutine ecat_prot(ecation_prot)
+! use calc_data
+! use comm_momo
+ integer i,j,k,subchap,itmp,inum
+ real(kind=8) :: xi,yi,zi,xj,yj,zj,ract,rcat0,epscalc,r06,r012,&
+ r7,r4,ecationcation
+ real(kind=8) xj_temp,yj_temp,zj_temp,xj_safe,yj_safe,zj_safe, &
+ dist_init,dist_temp,ecation_prot,rcal,rocal, &
+ Evan1,Evan2,EC,cm1mag,DASGL,delta,r0p,Epepcat, &
+ catl,cml,calpl, Etotal_p, Etotal_m,rtab,wdip,wmodquad,wquad1, &
+ wquad2,wvan1,E1,E2,wconst,wvan2,rcpm,dcmag,sin2thet,sinthet, &
+ costhet,v1m,v2m,wh2o,wc,rsecp,Ir,Irsecp,Irthrp,Irfourp,Irfiftp,&
+ Irsistp,Irseven,Irtwelv,Irthir,dE1dr,dE2dr,dEdcos,wquad2p,opt, &
+ rs,rthrp,rfourp,rsixp,reight,Irsixp,Ireight,Irtw,Irfourt, &
+ opt1,opt2,opt3,opt4,opt5,opt6,opt7,opt8,opt9,opt10,opt11,opt12,&
+ opt13,opt14,opt15,opt16,opt17,opt18,opt19, &
+ Equad1,Equad2,dscmag,v1dpv2,dscmag3,constA,constB,Edip,&
+ ndiv,ndivi
+ real(kind=8),dimension(3) ::dEvan1Cmcat,dEvan2Cmcat,dEeleccat,&
+ gg,r,EtotalCat,dEtotalCm,dEtotalCalp,dEvan1Cm,dEvan2Cm, &
+ dEtotalpep,dEtotalcat_num,dEddci,dEtotalcm_num,dEtotalcalp_num, &
+ tab1,tab2,tab3,diff,cm1,sc,p,tcat,talp,cm,drcp,drcp_norm,vcat, &
+ v1,v2,v3,myd_norm,dx,vcm,valpha,drdpep,dcosdpep,dcosddci,dEdpep,&
+ dEcCat,dEdipCm,dEdipCalp,dEquad1Cat,dEquad1Cm,dEquad1Calp, &
+ dEquad2Cat,dEquad2Cm,dEquad2Calpd,Evan1Cat,dEvan1Calp,dEvan2Cat,&
+ dEvan2Calp,dEtotalCat,dscvec,dEcCm,dEcCalp,dEdipCat,dEquad2Calp,&
+ dEvan1Cat
+ real(kind=8),dimension(6) :: vcatprm
+ ecation_prot=0.0d0
+! first lets calculate interaction with peptide groups
+ if (nres_molec(5).eq.0) return
+ itmp=0
+ do i=1,4
+ itmp=itmp+nres_molec(i)
+ enddo
+! do i=1,nres_molec(1)-1 ! loop over all peptide groups needs parralelization
+ do i=ibond_start,ibond_end
+! cycle
+ if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1)) cycle ! leave dummy atoms
+ xi=0.5d0*(c(1,i)+c(1,i+1))
+ yi=0.5d0*(c(2,i)+c(2,i+1))
+ zi=0.5d0*(c(3,i)+c(3,i+1))
+ xi=mod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=mod(yi,boxysize)
+ if (yi.lt.0) yi=yi+boxysize
+ zi=mod(zi,boxzsize)
+ if (zi.lt.0) zi=zi+boxzsize
+
+ do j=itmp+1,itmp+nres_molec(5)
+! print *,"WTF",itmp,j,i
+! all parameters were for Ca2+ to approximate single charge divide by two
+ ndiv=1.0
+ if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0
+ wconst=78*ndiv
+ wdip =1.092777950857032D2
+ wdip=wdip/wconst
+ wmodquad=-2.174122713004870D4
+ wmodquad=wmodquad/wconst
+ wquad1 = 3.901232068562804D1
+ wquad1=wquad1/wconst
+ wquad2 = 3
+ wquad2=wquad2/wconst
+ wvan1 = 0.1
+ wvan2 = 6
+! itmp=0
+
+ xj=c(1,j)
+ yj=c(2,j)
+ zj=c(3,j)
+ xj=dmod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=dmod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=dmod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
+! enddo
+! enddo
+ rcpm = sqrt(xj**2+yj**2+zj**2)
+ drcp_norm(1)=xj/rcpm
+ drcp_norm(2)=yj/rcpm
+ drcp_norm(3)=zj/rcpm
+ dcmag=0.0
+ do k=1,3
+ dcmag=dcmag+dc(k,i)**2
+ enddo
+ dcmag=dsqrt(dcmag)
do k=1,3
myd_norm(k)=dc(k,i)/dcmag
enddo
E2 = -wquad1*Irthrp*wquad2+wvan1*(wvan2**12*Irtwelv-&
2*wvan2**6*Irsistp)
ecation_prot = ecation_prot+E1+E2
+! print *,"ecatprot",i,j,ecation_prot,rcpm
dE1dr = -2*costhet*wdip*Irthrp-&
(4*wmodquad*Irfiftp+3*wquad1*Irfourp)*sin2thet
dE2dr = 3*wquad1*wquad2*Irfourp- &
enddo
cm1mag=sqrt(cm1(1)**2+cm1(2)**2+cm1(3)**2)
do j=itmp+1,itmp+nres_molec(5)
+ ndiv=1.0
+ if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0
+
xj=c(1,j)
yj=c(2,j)
zj=c(3,j)
endif
! enddo
! enddo
- if(itype(i,1).eq.15.or.itype(i,1).eq.16) then
+! 15- Glu 16-Asp
+ if((itype(i,1).eq.15.or.itype(i,1).eq.16).or.&
+ ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.&
+ (itype(i,1).eq.25))) then
if(itype(i,1).eq.16) then
inum=1
else
vcatprm(k)=catprm(k,inum)
enddo
dASGL=catprm(7,inum)
- do k=1,3
- vcm(k)=(cm1(k)/cm1mag)*dASGL+c(k,i+nres)
- valpha(k)=c(k,i)
- vcat(k)=c(k,j)
- enddo
- do k=1,3
+! do k=1,3
+! vcm(k)=(cm1(k)/cm1mag)*dASGL+c(k,i+nres)
+ vcm(1)=(cm1(1)/cm1mag)*dASGL+xi
+ vcm(2)=(cm1(2)/cm1mag)*dASGL+yi
+ vcm(3)=(cm1(3)/cm1mag)*dASGL+zi
+
+! valpha(k)=c(k,i)
+! vcat(k)=c(k,j)
+ if (subchap.eq.1) then
+ vcat(1)=xj_temp
+ vcat(2)=yj_temp
+ vcat(3)=zj_temp
+ else
+ vcat(1)=xj_safe
+ vcat(2)=yj_safe
+ vcat(3)=zj_safe
+ endif
+ valpha(1)=xi-c(1,i+nres)+c(1,i)
+ valpha(2)=yi-c(2,i+nres)+c(2,i)
+ valpha(3)=zi-c(3,i+nres)+c(3,i)
+
+! enddo
+ do k=1,3
dx(k) = vcat(k)-vcm(k)
enddo
do k=1,3
! The weights of the energy function calculated from
!The quantum mechanical GAMESS simulations of calcium with ASP/GLU
- wh2o=78
+ if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
+ ndivi=0.5
+ else
+ ndivi=1.0
+ endif
+ ndiv=1.0
+ if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0
+
+ wh2o=78*ndivi*ndiv
wc = vcatprm(1)
wc=wc/wh2o
wdip =vcatprm(2)
wquad1=wquad1/wh2o
wquad2 = vcatprm(4)
wquad2=wquad2/wh2o
- wquad2p = 1-wquad2
+ wquad2p = 1.0d0-wquad2
wvan1 = vcatprm(5)
wvan2 =vcatprm(6)
opt = dx(1)**2+dx(2)**2
rsixp = rfourp*rsecp
reight=rsixp*rsecp
Ir = 1.0d0/rs
- Irsecp = 1/rsecp
+ Irsecp = 1.0d0/rsecp
Irthrp = Irsecp/rs
Irfourp = Irthrp/rs
- Irsixp = 1/rsixp
- Ireight=1/reight
+ Irsixp = 1.0d0/rsixp
+ Ireight=1.0d0/reight
Irtw=Irsixp*Irsixp
Irthir=Irtw/rs
Irfourt=Irthir/rs
vcatprm(k)=catprm(k,inum)
enddo
dASGL=catprm(7,inum)
- do k=1,3
- vcm(k)=(cm1(k)/cm1mag)*dASGL+c(k,i+nres)
- valpha(k)=c(k,i)
- vcat(k)=c(k,j)
- enddo
+! do k=1,3
+! vcm(k)=(cm1(k)/cm1mag)*dASGL+c(k,i+nres)
+! valpha(k)=c(k,i)
+! vcat(k)=c(k,j)
+! enddo
+ vcm(1)=(cm1(1)/cm1mag)*dASGL+xi
+ vcm(2)=(cm1(2)/cm1mag)*dASGL+yi
+ vcm(3)=(cm1(3)/cm1mag)*dASGL+zi
+ if (subchap.eq.1) then
+ vcat(1)=xj_temp
+ vcat(2)=yj_temp
+ vcat(3)=zj_temp
+ else
+ vcat(1)=xj_safe
+ vcat(2)=yj_safe
+ vcat(3)=zj_safe
+ endif
+ valpha(1)=xi-c(1,i+nres)+c(1,i)
+ valpha(2)=yi-c(2,i+nres)+c(2,i)
+ valpha(3)=zi-c(3,i+nres)+c(3,i)
+
do k=1,3
dx(k) = vcat(k)-vcm(k)
v1dpv2 = v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3)
! The weights of the energy function calculated from
!The quantum mechanical GAMESS simulations of ASN/GLN with calcium
- wh2o=78
+ ndiv=1.0
+ if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0
+
+ wh2o=78*ndiv
wdip =vcatprm(2)
wdip=wdip/wh2o
wquad1 =vcatprm(3)
dscmag = 0.0d0
do k=1,3
dscvec(k) = c(k,i+nres)-c(k,i)
+! TU SPRAWDZ???
+! dscvec(1) = xj
+! dscvec(2) = yj
+! dscvec(3) = zj
+
dscmag = dscmag+dscvec(k)*dscvec(k)
enddo
dscmag3 = dscmag
else
rcal = 0.0d0
do k=1,3
- r(k) = c(k,j)-c(k,i+nres)
+! r(k) = c(k,j)-c(k,i+nres)
+ r(1) = xj
+ r(2) = yj
+ r(3) = zj
rcal = rcal+r(k)*r(k)
enddo
ract=sqrt(rcal)
!c! Compute head-head and head-tail energies for each state
isel = iabs(Qi) + iabs(Qj)
+! double charge for Phophorylated! itype - 25,27,27
+! if ((itype(i).eq.27).or.(itype(i).eq.26).or.(itype(i).eq.25)) then
+! Qi=Qi*2
+! Qij=Qij*2
+! endif
+! if ((itype(j).eq.27).or.(itype(j).eq.26).or.(itype(j).eq.25)) then
+! Qj=Qj*2
+! Qij=Qij*2
+! endif
+
! isel=0
IF (isel.eq.0) THEN
!c! No charges - do nothing
ELSE IF (isel.eq.1 .and. iabs(Qi).eq.1) THEN
!c! Charge-nonpolar interactions
+ if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
+ Qi=Qi*2
+ Qij=Qij*2
+ endif
+ if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then
+ Qj=Qj*2
+ Qij=Qij*2
+ endif
+
CALL eqn(epol)
eheadtail = epol
! eheadtail = 0.0d0
ELSE IF (isel.eq.1 .and. iabs(Qj).eq.1) THEN
!c! Nonpolar-charge interactions
+ if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
+ Qi=Qi*2
+ Qij=Qij*2
+ endif
+ if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then
+ Qj=Qj*2
+ Qij=Qij*2
+ endif
+
CALL enq(epol)
eheadtail = epol
! eheadtail = 0.0d0
ELSE IF (isel.eq.3 .and. icharge(itypj).eq.2) THEN
!c! Charge-dipole interactions
+ if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
+ Qi=Qi*2
+ Qij=Qij*2
+ endif
+ if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then
+ Qj=Qj*2
+ Qij=Qij*2
+ endif
+
CALL eqd(ecl, elj, epol)
eheadtail = ECL + elj + epol
! eheadtail = 0.0d0
ELSE IF (isel.eq.3 .and. icharge(itypi).eq.2) THEN
!c! Dipole-charge interactions
+ if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
+ Qi=Qi*2
+ Qij=Qij*2
+ endif
+ if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then
+ Qj=Qj*2
+ Qij=Qij*2
+ endif
CALL edq(ecl, elj, epol)
eheadtail = ECL + elj + epol
! eheadtail = 0.0d0
iabs(Qi).eq.1).and. &
nstate(itypi,itypj).eq.1) THEN
!c! Same charge-charge interaction ( +/+ or -/- )
+ if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
+ Qi=Qi*2
+ Qij=Qij*2
+ endif
+ if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then
+ Qj=Qj*2
+ Qij=Qij*2
+ endif
+
CALL eqq(Ecl,Egb,Epol,Fisocav,Elj)
eheadtail = ECL + Egb + Epol + Fisocav + Elj
! eheadtail = 0.0d0
iabs(Qi).eq.1).and. &
nstate(itypi,itypj).ne.1) THEN
!c! Different charge-charge interaction ( +/- or -/+ )
+ if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
+ Qi=Qi*2
+ Qij=Qij*2
+ endif
+ if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then
+ Qj=Qj*2
+ Qij=Qij*2
+ endif
+
CALL energy_quad(istate,eheadtail,Ecl,Egb,Epol,Fisocav,Elj,Equad)
END IF
END IF ! this endif ends the "catch the gly-gly" at the beggining of Fcav
use calc_data
use comm_momo
real (kind=8) :: facd3, facd4, federmaus, adler,&
- Ecl,Egb,Epol,Fisocav,Elj,Fgb
+ Ecl,Egb,Epol,Fisocav,Elj,Fgb,debkap
! integer :: k
!c! Epol and Gpol analytical parameters
alphapol1 = alphapol(itypi,itypj)
dGCLdOM12 = 0.0d0
ee0 = dexp(-( Rhead_sq ) / (4.0d0 * a12sq))
Fgb = sqrt( ( Rhead_sq ) + a12sq * ee0)
- Egb = -(332.0d0 * Qij * eps_inout_fac) / Fgb
+ debkap=debaykap(itypi,itypj)
+ Egb = -(332.0d0 * Qij *&
+ (1.0/eps_in-dexp(-debkap*Fgb)/eps_out)) / Fgb
! print *,"EGB WTF",Qij,eps_inout_fac,Fgb,itypi,itypj,eps_in,eps_out
!c! Derivative of Egb is Ggb...
- dGGBdFGB = -(-332.0d0 * Qij * eps_inout_fac) / (Fgb * Fgb)
+ dGGBdFGB = -(-332.0d0 * Qij * &
+ (1.0/eps_in-dexp(-debkap*Fgb)/eps_out))/(Fgb*Fgb)&
+ -(332.0d0 * Qij *&
+ (dexp(-debkap*Fgb)*debkap/eps_out))/ Fgb
dFGBdR = ( Rhead * ( 2.0d0 - (0.5d0 * ee0) ) )/ ( 2.0d0 * Fgb )
dGGBdR = dGGBdFGB * dFGBdR
!c!-------------------------------------------------------------------
END DO
RETURN
END SUBROUTINE eqq
+
+ SUBROUTINE eqq_cat(Ecl,Egb,Epol,Fisocav,Elj)
+ use calc_data
+ use comm_momo
+ real (kind=8) :: facd3, facd4, federmaus, adler,&
+ Ecl,Egb,Epol,Fisocav,Elj,Fgb,debkap
+! integer :: k
+!c! Epol and Gpol analytical parameters
+ alphapol1 = alphapolcat(itypi,itypj)
+ alphapol2 = alphapolcat(itypj,itypi)
+!c! Fisocav and Gisocav analytical parameters
+ al1 = alphisocat(1,itypi,itypj)
+ al2 = alphisocat(2,itypi,itypj)
+ al3 = alphisocat(3,itypi,itypj)
+ al4 = alphisocat(4,itypi,itypj)
+ csig = (1.0d0 &
+ / dsqrt(sigiso1cat(itypi, itypj)**2.0d0 &
+ + sigiso2cat(itypi,itypj)**2.0d0))
+!c!
+ pis = sig0headcat(itypi,itypj)
+ eps_head = epsheadcat(itypi,itypj)
+ Rhead_sq = Rhead * Rhead
+!c! R1 - distance between head of ith side chain and tail of jth sidechain
+!c! R2 - distance between head of jth side chain and tail of ith sidechain
+ R1 = 0.0d0
+ R2 = 0.0d0
+ DO k = 1, 3
+!c! Calculate head-to-tail distances needed by Epol
+ R1=R1+(ctail(k,2)-chead(k,1))**2
+ R2=R2+(chead(k,2)-ctail(k,1))**2
+ END DO
+!c! Pitagoras
+ R1 = dsqrt(R1)
+ R2 = dsqrt(R2)
+
+!c! R1 = dsqrt((Rtail**2)+((dtail(1,itypi,itypj)
+!c! & +dhead(1,1,itypi,itypj))**2))
+!c! R2 = dsqrt((Rtail**2)+((dtail(2,itypi,itypj)
+!c! & +dhead(2,1,itypi,itypj))**2))
+
+!c!-------------------------------------------------------------------
+!c! Coulomb electrostatic interaction
+ Ecl = (332.0d0 * Qij) / Rhead
+!c! derivative of Ecl is Gcl...
+ dGCLdR = (-332.0d0 * Qij ) / Rhead_sq
+ dGCLdOM1 = 0.0d0
+ dGCLdOM2 = 0.0d0
+ dGCLdOM12 = 0.0d0
+ ee0 = dexp(-( Rhead_sq ) / (4.0d0 * a12sq))
+ Fgb = sqrt( ( Rhead_sq ) + a12sq * ee0)
+ debkap=debaykapcat(itypi,itypj)
+ Egb = -(332.0d0 * Qij *&
+ (1.0/eps_in-dexp(-debkap*Fgb)/eps_out)) / Fgb
+! print *,"EGB WTF",Qij,eps_inout_fac,Fgb,itypi,itypj,eps_in,eps_out
+!c! Derivative of Egb is Ggb...
+ dGGBdFGB = -(-332.0d0 * Qij * &
+ (1.0/eps_in-dexp(-debkap*Fgb)/eps_out))/(Fgb*Fgb)&
+ -(332.0d0 * Qij *&
+ (dexp(-debkap*Fgb)*debkap/eps_out))/ Fgb
+ dFGBdR = ( Rhead * ( 2.0d0 - (0.5d0 * ee0) ) )/ ( 2.0d0 * Fgb )
+ dGGBdR = dGGBdFGB * dFGBdR
+!c!-------------------------------------------------------------------
+!c! Fisocav - isotropic cavity creation term
+!c! or "how much energy it costs to put charged head in water"
+ pom = Rhead * csig
+ top = al1 * (dsqrt(pom) + al2 * pom - al3)
+ bot = (1.0d0 + al4 * pom**12.0d0)
+ botsq = bot * bot
+ FisoCav = top / bot
+! write (*,*) "Rhead = ",Rhead
+! write (*,*) "csig = ",csig
+! write (*,*) "pom = ",pom
+! write (*,*) "al1 = ",al1
+! write (*,*) "al2 = ",al2
+! write (*,*) "al3 = ",al3
+! write (*,*) "al4 = ",al4
+! write (*,*) "top = ",top
+! write (*,*) "bot = ",bot
+!c! Derivative of Fisocav is GCV...
+ dtop = al1 * ((1.0d0 / (2.0d0 * dsqrt(pom))) + al2)
+ dbot = 12.0d0 * al4 * pom ** 11.0d0
+ dGCVdR = ((dtop * bot - top * dbot) / botsq) * csig
+!c!-------------------------------------------------------------------
+!c! Epol
+!c! Polarization energy - charged heads polarize hydrophobic "neck"
+ MomoFac1 = (1.0d0 - chi1 * sqom2)
+ MomoFac2 = (1.0d0 - chi2 * sqom1)
+ RR1 = ( R1 * R1 ) / MomoFac1
+ RR2 = ( R2 * R2 ) / MomoFac2
+ ee1 = exp(-( RR1 / (4.0d0 * a12sq) ))
+ ee2 = exp(-( RR2 / (4.0d0 * a12sq) ))
+ fgb1 = sqrt( RR1 + a12sq * ee1 )
+ fgb2 = sqrt( RR2 + a12sq * ee2 )
+ epol = 332.0d0 * eps_inout_fac * ( &
+ (( alphapol1 / fgb1 )**4.0d0)+((alphapol2/fgb2) ** 4.0d0 ))
+!c! epol = 0.0d0
+ dPOLdFGB1 = -(1328.0d0 * eps_inout_fac * alphapol1 ** 4.0d0)&
+ / (fgb1 ** 5.0d0)
+ dPOLdFGB2 = -(1328.0d0 * eps_inout_fac * alphapol2 ** 4.0d0)&
+ / (fgb2 ** 5.0d0)
+ dFGBdR1 = ( (R1 / MomoFac1)* ( 2.0d0 - (0.5d0 * ee1) ) )&
+ / ( 2.0d0 * fgb1 )
+ dFGBdR2 = ( (R2 / MomoFac2)* ( 2.0d0 - (0.5d0 * ee2) ) )&
+ / ( 2.0d0 * fgb2 )
+ dFGBdOM2 = (((R1 * R1 * chi1 * om2) / (MomoFac1 * MomoFac1))&
+ * ( 2.0d0 - 0.5d0 * ee1) ) / ( 2.0d0 * fgb1 )
+ dFGBdOM1 = (((R2 * R2 * chi2 * om1) / (MomoFac2 * MomoFac2))&
+ * ( 2.0d0 - 0.5d0 * ee2) ) / ( 2.0d0 * fgb2 )
+ dPOLdR1 = dPOLdFGB1 * dFGBdR1
+!c! dPOLdR1 = 0.0d0
+ dPOLdR2 = dPOLdFGB2 * dFGBdR2
+!c! dPOLdR2 = 0.0d0
+ dPOLdOM1 = dPOLdFGB2 * dFGBdOM1
+!c! dPOLdOM1 = 0.0d0
+ dPOLdOM2 = dPOLdFGB1 * dFGBdOM2
+!c! dPOLdOM2 = 0.0d0
+!c!-------------------------------------------------------------------
+!c! Elj
+!c! Lennard-Jones 6-12 interaction between heads
+ pom = (pis / Rhead)**6.0d0
+ Elj = 4.0d0 * eps_head * pom * (pom-1.0d0)
+!c! derivative of Elj is Glj
+ dGLJdR = 4.0d0 * eps_head*(((-12.0d0*pis**12.0d0)/(Rhead**13.0d0))&
+ + (( 6.0d0*pis**6.0d0) /(Rhead**7.0d0)))
+!c!-------------------------------------------------------------------
+!c! Return the results
+!c! These things do the dRdX derivatives, that is
+!c! allow us to change what we see from function that changes with
+!c! distance to function that changes with LOCATION (of the interaction
+!c! site)
+ DO k = 1, 3
+ erhead(k) = Rhead_distance(k)/Rhead
+ erhead_tail(k,1) = ((ctail(k,2)-chead(k,1))/R1)
+ erhead_tail(k,2) = ((chead(k,2)-ctail(k,1))/R2)
+ END DO
+
+ erdxi = scalar( erhead(1), dC_norm(1,i+nres) )
+ erdxj = scalar( erhead(1), dC_norm(1,j) )
+ bat = scalar( erhead_tail(1,1), dC_norm(1,i+nres) )
+ federmaus = scalar(erhead_tail(1,1),dC_norm(1,j))
+ eagle = scalar( erhead_tail(1,2), dC_norm(1,j) )
+ adler = scalar( erhead_tail(1,2), dC_norm(1,i+nres) )
+ facd1 = d1 * vbld_inv(i+nres)
+ facd2 = d2 * vbld_inv(j)
+ facd3 = dtailcat(1,itypi,itypj) * vbld_inv(i+nres)
+ facd4 = dtailcat(2,itypi,itypj) * vbld_inv(j)
+
+!c! Now we add appropriate partial derivatives (one in each dimension)
+ DO k = 1, 3
+ hawk = (erhead_tail(k,1) + &
+ facd1 * (erhead_tail(k,1) - bat * dC_norm(k,i+nres)))
+ condor = (erhead_tail(k,2) + &
+ facd2 * (erhead_tail(k,2) - eagle * dC_norm(k,j)))
+
+ pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
+ gvdwx(k,i) = gvdwx(k,i) &
+ - dGCLdR * pom&
+ - dGGBdR * pom&
+ - dGCVdR * pom&
+ - dPOLdR1 * hawk&
+ - dPOLdR2 * (erhead_tail(k,2)&
+ -facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres)))&
+ - dGLJdR * pom
+
+ pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j))
+ gvdwx(k,j) = gvdwx(k,j)+ dGCLdR * pom&
+ + dGGBdR * pom+ dGCVdR * pom&
+ + dPOLdR1 * (erhead_tail(k,1)&
+ -facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j)))&
+ + dPOLdR2 * condor + dGLJdR * pom
+
+ gvdwc(k,i) = gvdwc(k,i) &
+ - dGCLdR * erhead(k)&
+ - dGGBdR * erhead(k)&
+ - dGCVdR * erhead(k)&
+ - dPOLdR1 * erhead_tail(k,1)&
+ - dPOLdR2 * erhead_tail(k,2)&
+ - dGLJdR * erhead(k)
+
+ gvdwc(k,j) = gvdwc(k,j) &
+ + dGCLdR * erhead(k) &
+ + dGGBdR * erhead(k) &
+ + dGCVdR * erhead(k) &
+ + dPOLdR1 * erhead_tail(k,1) &
+ + dPOLdR2 * erhead_tail(k,2)&
+ + dGLJdR * erhead(k)
+
+ END DO
+ RETURN
+ END SUBROUTINE eqq_cat
!c!-------------------------------------------------------------------
SUBROUTINE energy_quad(istate,eheadtail,Ecl,Egb,Epol,Fisocav,Elj,Equad)
use comm_momo
facd4 = dtail(2,itypi,itypj) * vbld_inv(j+nres)
DO k = 1, 3
- hawk = (erhead_tail(k,1) + &
- facd1 * (erhead_tail(k,1) - bat * dC_norm(k,i+nres)))
+ hawk = (erhead_tail(k,1) + &
+ facd1 * (erhead_tail(k,1) - bat * dC_norm(k,i+nres)))
+
+ gvdwx(k,i) = gvdwx(k,i) &
+ - dPOLdR1 * hawk
+ gvdwx(k,j) = gvdwx(k,j) &
+ + dPOLdR1 * (erhead_tail(k,1) &
+ -facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j+nres)))
+
+ gvdwc(k,i) = gvdwc(k,i) - dPOLdR1 * erhead_tail(k,1)
+ gvdwc(k,j) = gvdwc(k,j) + dPOLdR1 * erhead_tail(k,1)
+
+ END DO
+ RETURN
+ END SUBROUTINE eqn
+ SUBROUTINE enq(Epol)
+ use calc_data
+ use comm_momo
+ double precision facd3, adler,epol
+ alphapol2 = alphapol(itypj,itypi)
+!c! R2 - distance between head of jth side chain and tail of ith sidechain
+ R2 = 0.0d0
+ DO k = 1, 3
+!c! Calculate head-to-tail distances
+ R2=R2+(chead(k,2)-ctail(k,1))**2
+ END DO
+!c! Pitagoras
+ R2 = dsqrt(R2)
+
+!c! R1 = dsqrt((Rtail**2)+((dtail(1,itypi,itypj)
+!c! & +dhead(1,1,itypi,itypj))**2))
+!c! R2 = dsqrt((Rtail**2)+((dtail(2,itypi,itypj)
+!c! & +dhead(2,1,itypi,itypj))**2))
+!c------------------------------------------------------------------------
+!c Polarization energy
+ MomoFac2 = (1.0d0 - chi2 * sqom1)
+ RR2 = R2 * R2 / MomoFac2
+ ee2 = exp(-(RR2 / (4.0d0 * a12sq)))
+ fgb2 = sqrt(RR2 + a12sq * ee2)
+ epol = 332.0d0 * eps_inout_fac * ((alphapol2/fgb2) ** 4.0d0 )
+ dPOLdFGB2 = -(1328.0d0 * eps_inout_fac * alphapol2 ** 4.0d0) &
+ / (fgb2 ** 5.0d0)
+ dFGBdR2 = ( (R2 / MomoFac2) &
+ * ( 2.0d0 - (0.5d0 * ee2) ) ) &
+ / (2.0d0 * fgb2)
+ dFGBdOM1 = (((R2 * R2 * chi2 * om1) / (MomoFac2 * MomoFac2)) &
+ * (2.0d0 - 0.5d0 * ee2) ) &
+ / (2.0d0 * fgb2)
+ dPOLdR2 = dPOLdFGB2 * dFGBdR2
+!c! dPOLdR2 = 0.0d0
+ dPOLdOM1 = dPOLdFGB2 * dFGBdOM1
+!c! dPOLdOM1 = 0.0d0
+ dPOLdOM2 = 0.0d0
+!c!-------------------------------------------------------------------
+!c! Return the results
+!c! (See comments in Eqq)
+ DO k = 1, 3
+ erhead_tail(k,2) = ((chead(k,2)-ctail(k,1))/R2)
+ END DO
+ eagle = scalar( erhead_tail(1,2), dC_norm(1,j+nres) )
+ adler = scalar( erhead_tail(1,2), dC_norm(1,i+nres) )
+ facd2 = d2 * vbld_inv(j+nres)
+ facd3 = dtail(1,itypi,itypj) * vbld_inv(i+nres)
+ DO k = 1, 3
+ condor = (erhead_tail(k,2) &
+ + facd2 * (erhead_tail(k,2) - eagle * dC_norm(k,j+nres)))
gvdwx(k,i) = gvdwx(k,i) &
- - dPOLdR1 * hawk
- gvdwx(k,j) = gvdwx(k,j) &
- + dPOLdR1 * (erhead_tail(k,1) &
- -facd4 * (erhead_tail(k,1) - federmaus * dC_norm(k,j+nres)))
+ - dPOLdR2 * (erhead_tail(k,2) &
+ -facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres)))
+ gvdwx(k,j) = gvdwx(k,j) &
+ + dPOLdR2 * condor
- gvdwc(k,i) = gvdwc(k,i) - dPOLdR1 * erhead_tail(k,1)
- gvdwc(k,j) = gvdwc(k,j) + dPOLdR1 * erhead_tail(k,1)
+ gvdwc(k,i) = gvdwc(k,i) &
+ - dPOLdR2 * erhead_tail(k,2)
+ gvdwc(k,j) = gvdwc(k,j) &
+ + dPOLdR2 * erhead_tail(k,2)
END DO
- RETURN
- END SUBROUTINE eqn
- SUBROUTINE enq(Epol)
+ RETURN
+ END SUBROUTINE enq
+
+ SUBROUTINE enq_cat(Epol)
use calc_data
use comm_momo
double precision facd3, adler,epol
- alphapol2 = alphapol(itypj,itypi)
+ alphapol2 = alphapolcat(itypj,itypi)
!c! R2 - distance between head of jth side chain and tail of ith sidechain
R2 = 0.0d0
DO k = 1, 3
dPOLdOM1 = dPOLdFGB2 * dFGBdOM1
!c! dPOLdOM1 = 0.0d0
dPOLdOM2 = 0.0d0
+
!c!-------------------------------------------------------------------
!c! Return the results
!c! (See comments in Eqq)
DO k = 1, 3
erhead_tail(k,2) = ((chead(k,2)-ctail(k,1))/R2)
END DO
- eagle = scalar( erhead_tail(1,2), dC_norm(1,j+nres) )
+ eagle = scalar( erhead_tail(1,2), dC_norm(1,j) )
adler = scalar( erhead_tail(1,2), dC_norm(1,i+nres) )
facd2 = d2 * vbld_inv(j+nres)
- facd3 = dtail(1,itypi,itypj) * vbld_inv(i+nres)
+ facd3 = dtailcat(1,itypi,itypj) * vbld_inv(i+nres)
DO k = 1, 3
condor = (erhead_tail(k,2) &
- + facd2 * (erhead_tail(k,2) - eagle * dC_norm(k,j+nres)))
+ + facd2 * (erhead_tail(k,2) - eagle * dC_norm(k,j)))
gvdwx(k,i) = gvdwx(k,i) &
- dPOLdR2 * (erhead_tail(k,2) &
END DO
RETURN
- END SUBROUTINE enq
+ END SUBROUTINE enq_cat
+
SUBROUTINE eqd(Ecl,Elj,Epol)
use calc_data
use comm_momo
END DO
RETURN
END SUBROUTINE edq
+
+ SUBROUTINE edq_cat(Ecl,Elj,Epol)
+ use comm_momo
+ use calc_data
+
+ double precision facd3, adler,ecl,elj,epol
+ alphapol2 = alphapolcat(itypj,itypi)
+ w1 = wqdipcat(1,itypi,itypj)
+ w2 = wqdipcat(2,itypi,itypj)
+ pis = sig0headcat(itypi,itypj)
+ eps_head = epsheadcat(itypi,itypj)
+!c!-------------------------------------------------------------------
+!c! R2 - distance between head of jth side chain and tail of ith sidechain
+ R2 = 0.0d0
+ DO k = 1, 3
+!c! Calculate head-to-tail distances
+ R2=R2+(chead(k,2)-ctail(k,1))**2
+ END DO
+!c! Pitagoras
+ R2 = dsqrt(R2)
+
+!c! R1 = dsqrt((Rtail**2)+((dtail(1,itypi,itypj)
+!c! & +dhead(1,1,itypi,itypj))**2))
+!c! R2 = dsqrt((Rtail**2)+((dtail(2,itypi,itypj)
+!c! & +dhead(2,1,itypi,itypj))**2))
+
+
+!c!-------------------------------------------------------------------
+!c! ecl
+ sparrow = w1 * Qi * om1
+ hawk = w2 * Qi * Qi * (1.0d0 - sqom2)
+ ECL = sparrow / Rhead**2.0d0 &
+ - hawk / Rhead**4.0d0
+!c!-------------------------------------------------------------------
+!c! derivative of ecl is Gcl
+!c! dF/dr part
+ dGCLdR = - 2.0d0 * sparrow / Rhead**3.0d0 &
+ + 4.0d0 * hawk / Rhead**5.0d0
+!c! dF/dom1
+ dGCLdOM1 = (w1 * Qi) / (Rhead**2.0d0)
+!c! dF/dom2
+ dGCLdOM2 = (2.0d0 * w2 * Qi * Qi * om2) / (Rhead ** 4.0d0)
+!c--------------------------------------------------------------------
+!c--------------------------------------------------------------------
+!c Polarization energy
+!c Epol
+ MomoFac2 = (1.0d0 - chi2 * sqom1)
+ RR2 = R2 * R2 / MomoFac2
+ ee2 = exp(-(RR2 / (4.0d0 * a12sq)))
+ fgb2 = sqrt(RR2 + a12sq * ee2)
+ epol = 332.0d0 * eps_inout_fac * ((alphapol2/fgb2) ** 4.0d0 )
+ dPOLdFGB2 = -(1328.0d0 * eps_inout_fac * alphapol2 ** 4.0d0) &
+ / (fgb2 ** 5.0d0)
+ dFGBdR2 = ( (R2 / MomoFac2) &
+ * ( 2.0d0 - (0.5d0 * ee2) ) ) &
+ / (2.0d0 * fgb2)
+ dFGBdOM1 = (((R2 * R2 * chi2 * om1) / (MomoFac2 * MomoFac2)) &
+ * (2.0d0 - 0.5d0 * ee2) ) &
+ / (2.0d0 * fgb2)
+ dPOLdR2 = dPOLdFGB2 * dFGBdR2
+!c! dPOLdR2 = 0.0d0
+ dPOLdOM1 = dPOLdFGB2 * dFGBdOM1
+!c! dPOLdOM1 = 0.0d0
+ dPOLdOM2 = 0.0d0
+!c!-------------------------------------------------------------------
+!c! Elj
+ pom = (pis / Rhead)**6.0d0
+ Elj = 4.0d0 * eps_head * pom * (pom-1.0d0)
+!c! derivative of Elj is Glj
+ dGLJdR = 4.0d0 * eps_head &
+ * (((-12.0d0*pis**12.0d0)/(Rhead**13.0d0)) &
+ + (( 6.0d0*pis**6.0d0) /(Rhead**7.0d0)))
+!c!-------------------------------------------------------------------
+
+!c! Return the results
+!c! (see comments in Eqq)
+ DO k = 1, 3
+ erhead(k) = Rhead_distance(k)/Rhead
+ erhead_tail(k,2) = ((chead(k,2)-ctail(k,1))/R2)
+ END DO
+ erdxi = scalar( erhead(1), dC_norm(1,i+nres) )
+ erdxj = scalar( erhead(1), dC_norm(1,j) )
+ eagle = scalar( erhead_tail(1,2), dC_norm(1,j) )
+ adler = scalar( erhead_tail(1,2), dC_norm(1,i+nres) )
+ facd1 = d1 * vbld_inv(i+nres)
+ facd2 = d2 * vbld_inv(j)
+ facd3 = dtail(1,itypi,itypj) * vbld_inv(i+nres)
+ DO k = 1, 3
+ condor = (erhead_tail(k,2) &
+ + facd2 * (erhead_tail(k,2) - eagle * dC_norm(k,j)))
+
+ pom = erhead(k)+facd1*(erhead(k)-erdxi*dC_norm(k,i+nres))
+ gvdwx(k,i) = gvdwx(k,i) &
+ - dGCLdR * pom &
+ - dPOLdR2 * (erhead_tail(k,2) &
+ -facd3 * (erhead_tail(k,2) - adler * dC_norm(k,i+nres))) &
+ - dGLJdR * pom
+
+ pom = erhead(k)+facd2*(erhead(k)-erdxj*dC_norm(k,j))
+ gvdwx(k,j) = gvdwx(k,j) &
+ + dGCLdR * pom &
+ + dPOLdR2 * condor &
+ + dGLJdR * pom
+
+
+ gvdwc(k,i) = gvdwc(k,i) &
+ - dGCLdR * erhead(k) &
+ - dPOLdR2 * erhead_tail(k,2) &
+ - dGLJdR * erhead(k)
+
+ gvdwc(k,j) = gvdwc(k,j) &
+ + dGCLdR * erhead(k) &
+ + dPOLdR2 * erhead_tail(k,2) &
+ + dGLJdR * erhead(k)
+
+ END DO
+ RETURN
+ END SUBROUTINE edq_cat
+
+
SUBROUTINE edd(ECL)
! IMPLICIT NONE
use comm_momo
dPOLdOM2 = 0.0d0
RETURN
END SUBROUTINE elgrad_init
+
+
+ SUBROUTINE elgrad_init_cat(eheadtail,Egb,Ecl,Elj,Equad,Epol)
+ use comm_momo
+ use calc_data
+ real(kind=8) :: eheadtail,Egb,Ecl,Elj,Equad,Epol,Rb
+ eps_out=80.0d0
+ itypi = itype(i,1)
+ itypj = itype(j,1)
+!c! 1/(Gas Constant * Thermostate temperature) = BetaT
+!c! ENABLE THIS LINE WHEN USING CHECKGRAD!!!
+!c! t_bath = 300
+!c! BetaT = 1.0d0 / (t_bath * Rb)i
+ Rb=0.001986d0
+ BetaT = 1.0d0 / (298.0d0 * Rb)
+!c! Gay-berne var's
+ sig0ij = sigmacat( itypi,itypj )
+ chi1 = chicat( itypi, itypj )
+! chi2 = chi( itypj, itypi )
+ chi2 = 0.0d0
+! chi12 = chi1 * chi2
+ chi12 = 0.0d0
+ chip1 = chippcat( itypi, itypj )
+! chip2 = chipp( itypj, itypi )
+ chip2 = 0.0d0
+! chip12 = chip1 * chip2
+ chip12 = 0.0d0
+! chi1=0.0
+! chi2=0.0
+! chi12=0.0
+! chip1=0.0
+! chip2=0.0
+! chip12=0.0
+!c! not used by momo potential, but needed by sc_angular which is shared
+!c! by all energy_potential subroutines
+ alf1 = 0.0d0
+ alf2 = 0.0d0
+ alf12 = 0.0d0
+!c! location, location, location
+! xj = c( 1, nres+j ) - xi
+! yj = c( 2, nres+j ) - yi
+! zj = c( 3, nres+j ) - zi
+ dxj = dc_norm( 1, nres+j )
+ dyj = dc_norm( 2, nres+j )
+ dzj = dc_norm( 3, nres+j )
+!c! distance from center of chain(?) to polar/charged head
+ d1 = dheadcat(1, 1, itypi, itypj)
+ d2 = dheadcat(2, 1, itypi, itypj)
+!c! ai*aj from Fgb
+ a12sq = rborncat(itypi,itypj) * rborncat(itypj,itypi)
+!c! a12sq = a12sq * a12sq
+!c! charge of amino acid itypi is...
+ Qi = ichargecat(itypi)
+ Qj = ichargecat(itypj)
+ Qij = Qi * Qj
+!c! chis1,2,12
+ chis1 = chiscat(itypi,itypj)
+! chis2 = chis(itypj,itypi)
+ chis2 = 0.0d0
+! chis12 = chis1 * chis2
+ chis12 = 0.0d0
+ sig1 = sigmap1cat(itypi,itypj)
+ sig2 = sigmap2cat(itypi,itypj)
+!c! alpha factors from Fcav/Gcav
+ b1cav = alphasurcat(1,itypi,itypj)
+! b1cav=0.0
+ b2cav = alphasurcat(2,itypi,itypj)
+ b3cav = alphasurcat(3,itypi,itypj)
+ b4cav = alphasurcat(4,itypi,itypj)
+ wqd = wquadcat(itypi, itypj)
+!c! used by Fgb
+ eps_in = epsintabcat(itypi,itypj)
+ eps_inout_fac = ( (1.0d0/eps_in) - (1.0d0/eps_out))
+!c!-------------------------------------------------------------------
+!c! tail location and distance calculations
+ Rtail = 0.0d0
+ DO k = 1, 3
+ ctail(k,1)=c(k,i+nres)-dtailcat(1,itypi,itypj)*dc_norm(k,nres+i)
+ ctail(k,2)=c(k,j+nres)-dtailcat(2,itypi,itypj)*dc_norm(k,nres+j)
+ END DO
+!c! tail distances will be themselves usefull elswhere
+!c1 (in Gcav, for example)
+ Rtail_distance(1) = ctail( 1, 2 ) - ctail( 1,1 )
+ Rtail_distance(2) = ctail( 2, 2 ) - ctail( 2,1 )
+ Rtail_distance(3) = ctail( 3, 2 ) - ctail( 3,1 )
+ Rtail = dsqrt( &
+ (Rtail_distance(1)*Rtail_distance(1)) &
+ + (Rtail_distance(2)*Rtail_distance(2)) &
+ + (Rtail_distance(3)*Rtail_distance(3)))
+!c!-------------------------------------------------------------------
+!c! Calculate location and distance between polar heads
+!c! distance between heads
+!c! for each one of our three dimensional space...
+ d1 = dheadcat(1, 1, itypi, itypj)
+ d2 = dheadcat(2, 1, itypi, itypj)
+
+ DO k = 1,3
+!c! location of polar head is computed by taking hydrophobic centre
+!c! and moving by a d1 * dc_norm vector
+!c! see unres publications for very informative images
+ chead(k,1) = c(k, i+nres) + d1 * dc_norm(k, i+nres)
+ chead(k,2) = c(k, j+nres) + d2 * dc_norm(k, j+nres)
+!c! distance
+!c! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres))
+!c! Rsc(k) = Rsc_distance(k) * Rsc_distance(k)
+ Rhead_distance(k) = chead(k,2) - chead(k,1)
+ END DO
+!c! pitagoras (root of sum of squares)
+ Rhead = dsqrt( &
+ (Rhead_distance(1)*Rhead_distance(1)) &
+ + (Rhead_distance(2)*Rhead_distance(2)) &
+ + (Rhead_distance(3)*Rhead_distance(3)))
+!c!-------------------------------------------------------------------
+!c! zero everything that should be zero'ed
+ Egb = 0.0d0
+ ECL = 0.0d0
+ Elj = 0.0d0
+ Equad = 0.0d0
+ Epol = 0.0d0
+ eheadtail = 0.0d0
+ dGCLdOM1 = 0.0d0
+ dGCLdOM2 = 0.0d0
+ dGCLdOM12 = 0.0d0
+ dPOLdOM1 = 0.0d0
+ dPOLdOM2 = 0.0d0
+ RETURN
+ END SUBROUTINE elgrad_init_cat
+
+
+ double precision function tschebyshev(m,n,x,y)
+ implicit none
+ integer i,m,n
+ double precision x(n),y,yy(0:maxvar),aux
+!c Tschebyshev polynomial. Note that the first term is omitted
+!c m=0: the constant term is included
+!c m=1: the constant term is not included
+ yy(0)=1.0d0
+ yy(1)=y
+ do i=2,n
+ yy(i)=2*yy(1)*yy(i-1)-yy(i-2)
+ enddo
+ aux=0.0d0
+ do i=m,n
+ aux=aux+x(i)*yy(i)
+ enddo
+ tschebyshev=aux
+ return
+ end function tschebyshev
+!C--------------------------------------------------------------------------
+ double precision function gradtschebyshev(m,n,x,y)
+ implicit none
+ integer i,m,n
+ double precision x(n+1),y,yy(0:maxvar),aux
+!c Tschebyshev polynomial. Note that the first term is omitted
+!c m=0: the constant term is included
+!c m=1: the constant term is not included
+ yy(0)=1.0d0
+ yy(1)=2.0d0*y
+ do i=2,n
+ yy(i)=2*y*yy(i-1)-yy(i-2)
+ enddo
+ aux=0.0d0
+ do i=m,n
+ aux=aux+x(i+1)*yy(i)*(i+1)
+!C print *, x(i+1),yy(i),i
+ enddo
+ gradtschebyshev=aux
+ return
+ end function gradtschebyshev
+
+
+
+
+
end module energy
projects / unres4.git / blobdiff