gvdwpp_nucl
!-----------------------------NUCLEIC-PROTEIN GRADIENT
real(kind=8),dimension(:,:),allocatable :: gvdwx_scbase,gvdwc_scbase,&
- gvdwx_pepbase,gvdwc_pepbase,gvdwx_scpho,gvdwc_scpho
+ gvdwx_pepbase,gvdwc_pepbase,gvdwx_scpho,gvdwc_scpho,&
+ gvdwc_peppho
!------------------------------IONS GRADIENT
real(kind=8),dimension(:,:),allocatable :: gradcatcat, &
gradpepcat,gradpepcatx
! energies for ions
real(kind=8) :: ecation_prot,ecationcation
! energies for protein nucleic acid interaction
- real(kind=8) :: escbase,epepbase,escpho
+ real(kind=8) :: escbase,epepbase,escpho,epeppho
#ifdef MPI
real(kind=8) :: weights_(n_ene) !,time_Bcast,time_Bcastw
weights_(42)=wcatprot
weights_(46)=wscbase
weights_(47)=wscpho
+ weights_(48)=wpeppho
! wcatcat= weights(41)
! wcatprot=weights(42)
wcatprot=weights(42)
wscbase=weights(46)
wscpho=weights(47)
+ wpeppho=weights(48)
endif
time_Bcast=time_Bcast+MPI_Wtime()-time00
time_Bcastw=time_Bcastw+MPI_Wtime()-time00
etube=0.0d0
endif
!--------------------------------------------------------
+! write (iout,*) "NRES_MOLEC(2),",nres_molec(2)
! print *,"before",ees,evdw1,ecorr
+ if (nres_molec(2).gt.0) then
call ebond_nucl(estr_nucl)
call ebend_nucl(ebe_nucl)
call etor_nucl(etors_nucl)
call epsb(evdwpsb,eelpsb)
call esb(esbloc)
call multibody_hb_nucl(ecorr_nucl,ecorr3_nucl,n_corr,n_corr1)
+ else
+ etors_nucl=0.0d0
+ estr_nucl=0.0d0
+ ebe_nucl=0.0d0
+ evdwsb=0.0d0
+ eelsb=0.0d0
+ esbloc=0.0d0
+ endif
+ if (nfgtasks.gt.1) then
+ if (fg_rank.eq.0) then
+ call ecatcat(ecationcation)
+ endif
+ else
call ecatcat(ecationcation)
+ endif
call ecat_prot(ecation_prot)
+ if (nres_molec(2).gt.0) then
call eprot_sc_base(escbase)
call epep_sc_base(epepbase)
call eprot_sc_phosphate(escpho)
+ call eprot_pep_phosphate(epeppho)
+ endif
! call ecatcat(ecationcation)
! print *,"after ebend", ebe_nucl
#ifdef TIMING
energia(46)=escbase
energia(47)=epepbase
energia(48)=escpho
+ energia(49)=epeppho
call sum_energy(energia,.true.)
if (dyn_ss) call dyn_set_nss
! print *," Processor",myrank," left SUM_ENERGY"
ebe_nucl,esbloc,etors_nucl,etors_d_nucl,ecorr_nucl,&
ecorr3_nucl
real(kind=8) :: ecation_prot,ecationcation
- real(kind=8) :: escbase,epepbase,escpho
+ real(kind=8) :: escbase,epepbase,escpho,epeppho
integer :: i
#ifdef MPI
integer :: ierr
escbase=energia(46)
epepbase=energia(47)
escpho=energia(48)
+ epeppho=energia(49)
! 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
+ +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho
#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
+ +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho
#endif
energia(0)=etot
! detecting NaNQ
ebe_nucl,esbloc,etors_nucl,etors_d_nucl,ecorr_nucl,&
ecorr3_nucl
real(kind=8) :: ecation_prot,ecationcation
- real(kind=8) :: escbase,epepbase,escpho
+ real(kind=8) :: escbase,epepbase,escpho,epeppho
etot=energia(0)
evdw=energia(1)
escbase=energia(46)
epepbase=energia(47)
escpho=energia(48)
+ epeppho=energia(49)
#ifdef SPLITELE
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,&
estr,wbond,ebe,wang,&
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,&
+ escbase,wscbase,epepbase,wpepbase,escpho,wscpho,epeppho,wpeppho,&
etot
10 format (/'Virtual-chain energies:'// &
'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ &
'ESCBASE=',1pE16.6,' WEIGHT=',1pD16.6,'(sc-prot nucl-base)'/ &
'EPEPBASE=',1pE16.6,' WEIGHT=',1pD16.6,'(pep-prot nucl-base)'/ &
'ESCPHO=',1pE16.6,' WEIGHT=',1pD16.6,'(sc-prot nucl-phosphate)'/&
+ 'EPEPPHO=',1pE16.6,' WEIGHT=',1pD16.6,'(pep-prot nucl-phosphate)'/&
'ETOT= ',1pE16.6,' (total)')
#else
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,&
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,&
+ escbase,wscbase,epepbase,wpepbase,escpho,wscpho,epeppho,wpeppho,&
etot
10 format (/'Virtual-chain energies:'// &
'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ &
'ESCBASE=',1pE16.6,' WEIGHT=',1pD16.6,'(sc-prot nucl-base)'/ &
'EPEPBASE=',1pE16.6,' WEIGHT=',1pD16.6,'(pep-prot nucl-base)'/ &
'ESCPHO=',1pE16.6,' WEIGHT=',1pD16.6,'(sc-prot nucl-phosphate)'/&
+ 'EPEPPHO=',1pE16.6,' WEIGHT=',1pD16.6,'(pep-prot nucl-phosphate)'/&
'ETOT= ',1pE16.6,' (total)')
#endif
return
endif
! if (i.gt. iatel_s+2 .and. i.lt.iatel_e+5) then
if (i.gt. nnt+2 .and. i.lt.nct+2) then
+ if (itype(i-2,1).eq.0) then
+ iti=ntortyp+1
+ else
iti = itortyp(itype(i-2,1))
+ endif
else
iti=ntortyp+1
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
+ else
iti1 = itortyp(itype(i-1,1))
+ endif
else
iti1=ntortyp+1
endif
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).le.ntyp) then
+ if (itype(i-1,1).eq.0) then
+ iti1=ntortyp+1
+ elseif (itype(i-1,1).le.ntyp) then
iti1 = itortyp(itype(i-1,1))
else
iti1=ntortyp+1
wcatcat*gradcatcat(j,i)+ &
wscbase*gvdwc_scbase(j,i)+ &
wpepbase*gvdwc_pepbase(j,i)+&
- wscpho*gvdwc_scpho(j,i)
+ wscpho*gvdwc_scpho(j,i)+ &
+ wpeppho*gvdwc_peppho(j,i)
+
wcatcat*gradcatcat(j,i)+ &
wscbase*gvdwc_scbase(j,i) &
wpepbase*gvdwc_pepbase(j,i)+&
- wscpho*gvdwc_scpho(j,i)
+ wscpho*gvdwc_scpho(j,i)+&
+ wpeppho*gvdwc_peppho(j,i)
enddo
+wscbase*gvdwx_scbase(j,i) &
+wpepbase*gvdwx_pepbase(j,i)&
+wscpho*gvdwx_scpho(j,i)
- if (i.eq.3) print *,"tu?", wscpho,gvdwx_scpho(j,i)
+! if (i.eq.3) print *,"tu?", wscpho,gvdwx_scpho(j,i)
enddo
enddo
gvdwc_pepbase(j,i)=0.0d0
gvdwx_scpho(j,i)=0.0d0
gvdwc_scpho(j,i)=0.0d0
+ gvdwc_peppho(j,i)=0.0d0
enddo
enddo
do i=0,nres
allocate(gvdwc_pepbase(3,-1:nres))
allocate(gvdwx_scpho(3,-1:nres))
allocate(gvdwc_scpho(3,-1:nres))
+ allocate(gvdwc_peppho(3,-1:nres))
allocate(dtheta(3,2,-1:nres))
!(3,2,maxres)
evdwij=evdwij*eps2rt
evdwsb=evdwsb+evdwij
if (lprn) then
- sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
- epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+ sigm=dabs(aa_nucl(itypi,itypj)/bb_nucl(itypi,itypj))**(1.0D0/6.0D0)
+ epsi=bb_nucl(itypi,itypj)**2/aa_nucl(itypi,itypj)
write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
restyp(itypi,2),i,restyp(itypj,2),j, &
epsi,sigm,chi1,chi2,chip1,chip2, &
r012 = r06**2
k0 = 332.0*(2.0*2.0)/80.0
itmp=0
+
do i=1,4
itmp=itmp+nres_molec(i)
enddo
- do i=itmp+1,itmp+nres_molec(i)-1
+ write(iout,*) "itmp",itmp
+ do i=itmp+1,itmp+nres_molec(5)-1
xi=c(1,i)
yi=c(2,i)
zi=c(3,i)
+
xi=mod(xi,boxxsize)
if (xi.lt.0) xi=xi+boxxsize
yi=mod(yi,boxysize)
if (yj.lt.0) yj=yj+boxysize
zj=dmod(zj,boxzsize)
if (zj.lt.0) zj=zj+boxzsize
+ write(iout,*) c(1,i),xi,xj,"xy",boxxsize
dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
xj_safe=xj
yj_safe=yj
gradcatcat(k,j)=gradcatcat(k,j)+gg(k)
enddo
+ write(iout,*) "ecatcat",i,j, ecationcation,xj,yj,zj
ecationcation=ecationcation+Evan1cat+Evan2cat+Eeleccat
enddo
enddo
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=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))
enddo ! j
enddo ! i
!------------------------------------------sidechains
- do i=1,nres_molec(1)
+! do i=1,nres_molec(1)
+ do i=ibond_start,ibond_end
if ((itype(i,1).eq.ntyp1)) cycle ! leave dummy atoms
! cycle
! print *,i,ecation_prot
real (kind=8),dimension(4):: ener
real(kind=8) :: b1,b2,b3,b4,egb,eps_in,eps_inout_fac,eps_out
real(kind=8) :: ECL,Elj,Equad,Epol,eheadtail,rhead,dGCLOM2,&
- sqom1,sqom2,sqom12,c1,c2,pom,Lambf,sparrow,&
+ sqom1,sqom2,sqom12,c1,c2,c3,pom,Lambf,sparrow,&
Chif,ChiLambf,bat,eagle,top,bot,botsq,Fcav,dtop,dFdR,dFdOM1,&
- dFdOM2,w1,w2,dGCLdR,dFdL,dFdOM12,dbot ,&
+ dFdOM2,w1,w2,w3,dGCLdR,dFdL,dFdOM12,dbot ,&
r1,eps_head,alphapol1,pis,facd2,d2,facd1,d1,erdxj,erdxi,federmaus,&
dPOLdR1,dFGBdOM2,dFGBdR1,dPOLdFGB1,RR1,MomoFac1,hawk,d1i,d1j,&
sig1,sig2,chis12,chis2,ee1,fgb1,a12sq,chis1
integer troll
eps_out=80.0d0
escbase=0.0d0
- do i=1,nres_molec(1)
+! do i=1,nres_molec(1)
+ do i=ibond_start,ibond_end
if (itype(i,1).eq.ntyp1_molec(1)) cycle
itypi = itype(i,1)
dxi = dc_norm(1,nres+i)
!Now dipole-dipole
if (wdipdip_scbase(2,itypi,itypj).gt.0.0d0) then
w1 = wdipdip_scbase(1,itypi,itypj)
- w2 = wdipdip_scbase(2,itypi,itypj)
+ w2 = -wdipdip_scbase(3,itypi,itypj)/2.0
+ w3 = wdipdip_scbase(2,itypi,itypj)
!c!-------------------------------------------------------------------
!c! ECL
fac = (om12 - 3.0d0 * om1 * om2)
c1 = (w1 / (Rhead**3.0d0)) * fac
c2 = (w2 / Rhead ** 6.0d0) &
* (4.0d0 + fac * fac -3.0d0 * (sqom1 + sqom2))
- ECL = c1 - c2
+ c3= (w3/ Rhead ** 6.0d0) &
+ * (2.0d0 - 2.0d0*fac*fac +3.0d0*(sqom1 + sqom2))
+ ECL = c1 - c2 + c3
!c! write (*,*) "w1 = ", w1
!c! write (*,*) "w2 = ", w2
!c! write (*,*) "om1 = ", om1
c1 = (-3.0d0 * w1 * fac) / (Rhead ** 4.0d0)
c2 = (-6.0d0 * w2) / (Rhead ** 7.0d0) &
* (4.0d0 + fac * fac - 3.0d0 * (sqom1 + sqom2))
- dGCLdR = c1 - c2
+ c3= (-6.0d0 * w3) / (Rhead ** 7.0d0) &
+ * (2.0d0 - 2.0d0*fac*fac +3.0d0*(sqom1 + sqom2))
+ dGCLdR = c1 - c2 + c3
!c! dECL/dom1
c1 = (-3.0d0 * w1 * om2 ) / (Rhead**3.0d0)
c2 = (-6.0d0 * w2) / (Rhead**6.0d0) &
* ( om2 * om12 - 3.0d0 * om1 * sqom2 + om1 )
- dGCLdOM1 = c1 - c2
+ c3 =(6.0d0*w3/ Rhead ** 6.0d0)*(om1-2.0d0*(fac)*(-om2))
+ dGCLdOM1 = c1 - c2 + c3
!c! dECL/dom2
c1 = (-3.0d0 * w1 * om1 ) / (Rhead**3.0d0)
c2 = (-6.0d0 * w2) / (Rhead**6.0d0) &
* ( om1 * om12 - 3.0d0 * sqom1 * om2 + om2 )
- dGCLdOM2 = c1 - c2
+ c3 =(6.0d0*w3/ Rhead ** 6.0d0)*(om2-2.0d0*(fac)*(-om1))
+ dGCLdOM2 = c1 - c2 + c3
!c! dECL/dom12
c1 = w1 / (Rhead ** 3.0d0)
c2 = ( 2.0d0 * w2 * fac ) / Rhead ** 6.0d0
- dGCLdOM12 = c1 - c2
+ c3 = (w3/ Rhead ** 6.0d0)*(-4.0d0*fac)
+ dGCLdOM12 = c1 - c2 + c3
DO k= 1, 3
erhead(k) = Rhead_distance(k)/Rhead
END DO
real (kind=8),dimension(4):: ener
real(kind=8) :: b1,b2,b3,b4,egb,eps_in,eps_inout_fac,eps_out
real(kind=8) :: ECL,Elj,Equad,Epol,eheadtail,rhead,dGCLOM2,&
- sqom1,sqom2,sqom12,c1,c2,pom,Lambf,sparrow,&
+ sqom1,sqom2,sqom12,c1,c2,c3,pom,Lambf,sparrow,&
Chif,ChiLambf,bat,eagle,top,bot,botsq,Fcav,dtop,dFdR,dFdOM1,&
- dFdOM2,w1,w2,dGCLdR,dFdL,dFdOM12,dbot ,&
+ dFdOM2,w1,w2,w3,dGCLdR,dFdL,dFdOM12,dbot ,&
r1,eps_head,alphapol1,pis,facd2,d2,facd1,d1,erdxj,erdxi,federmaus,&
dPOLdR1,dFGBdOM2,dFGBdR1,dPOLdFGB1,RR1,MomoFac1,hawk,d1i,d1j,&
sig1,sig2,chis12,chis2,ee1,fgb1,a12sq,chis1
integer troll
eps_out=80.0d0
epepbase=0.0d0
- do i=1,nres_molec(1)-1
+! do i=1,nres_molec(1)-1
+ do i=ibond_start,ibond_end
if (itype(i,1).eq.ntyp1_molec(1).or.itype(i+1,1).eq.ntyp1_molec(1)) cycle
!C itypi = itype(i,1)
dxi = dc_norm(1,i)
w1 = wdipdip_pepbase(1,itypj)
- w2 = wdipdip_pepbase(2,itypj)
+ w2 = -wdipdip_pepbase(3,itypj)/2.0
+ w3 = wdipdip_pepbase(2,itypj)
! w1=0.0d0
! w2=0.0d0
!c!-------------------------------------------------------------------
!c! ECL
+! w3=0.0d0
fac = (om12 - 3.0d0 * om1 * om2)
c1 = (w1 / (Rhead**3.0d0)) * fac
c2 = (w2 / Rhead ** 6.0d0) &
* (4.0d0 + fac * fac -3.0d0 * (sqom1 + sqom2))
- ECL = c1 - c2
+ c3= (w3/ Rhead ** 6.0d0) &
+ * (2.0d0 - 2.0d0*fac*fac +3.0d0*(sqom1 + sqom2))
+
+ ECL = c1 - c2 + c3
+
c1 = (-3.0d0 * w1 * fac) / (Rhead ** 4.0d0)
c2 = (-6.0d0 * w2) / (Rhead ** 7.0d0) &
* (4.0d0 + fac * fac - 3.0d0 * (sqom1 + sqom2))
- dGCLdR = c1 - c2
+ c3= (-6.0d0 * w3) / (Rhead ** 7.0d0) &
+ * (2.0d0 - 2.0d0*fac*fac +3.0d0*(sqom1 + sqom2))
+
+ dGCLdR = c1 - c2 + c3
!c! dECL/dom1
c1 = (-3.0d0 * w1 * om2 ) / (Rhead**3.0d0)
c2 = (-6.0d0 * w2) / (Rhead**6.0d0) &
* ( om2 * om12 - 3.0d0 * om1 * sqom2 + om1 )
- dGCLdOM1 = c1 - c2
+ c3 =(6.0d0*w3/ Rhead ** 6.0d0)*(om1-2.0d0*(fac)*(-om2))
+ dGCLdOM1 = c1 - c2 + c3
!c! dECL/dom2
c1 = (-3.0d0 * w1 * om1 ) / (Rhead**3.0d0)
c2 = (-6.0d0 * w2) / (Rhead**6.0d0) &
* ( om1 * om12 - 3.0d0 * sqom1 * om2 + om2 )
- dGCLdOM2 = c1 - c2
+ c3 =(6.0d0*w3/ Rhead ** 6.0d0)*(om2-2.0d0*(fac)*(-om1))
+
+ dGCLdOM2 = c1 - c2 + c3
!c! dECL/dom12
c1 = w1 / (Rhead ** 3.0d0)
c2 = ( 2.0d0 * w2 * fac ) / Rhead ** 6.0d0
- dGCLdOM12 = c1 - c2
+ c3 = (w3/ Rhead ** 6.0d0)*(-4.0d0*fac)
+ dGCLdOM12 = c1 - c2 + c3
DO k= 1, 3
erhead(k) = Rhead_distance(k)/Rhead
END DO
real(kind=8) :: evdw,sig0ij
real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
dist_temp, dist_init,aa,bb,ssgradlipi,ssgradlipj, &
- sslipi,sslipj,faclip
+ sslipi,sslipj,faclip,alpha_sco
integer :: ii
real(kind=8) :: fracinbuf
real (kind=8) :: escpho
integer troll
eps_out=80.0d0
escpho=0.0d0
- do i=1,nres_molec(1)
+! do i=1,nres_molec(1)
+ do i=ibond_start,ibond_end
if (itype(i,1).eq.ntyp1_molec(1)) cycle
itypi = itype(i,1)
dxi = dc_norm(1,nres+i)
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
dxj = dc_norm( 1,j )
dyj = dc_norm( 2,j )
dzj = dc_norm( 3,j )
+ dscj_inv = vbld_inv(j+1)
+
! Gay-berne var's
sig0ij = sigma_scpho(itypi )
chi1 = chi_scpho(itypi,1 )
!c! write (*,*) "Gvdwc(",k,",",i,")=", gvdwc(k,i)
!c! write (*,*) "Gvdwc(",k,",",j,")=", gvdwc(k,j)
! alphapol1 = alphapol_scpho(itypi)
- if (wqq_scpho(itypi).gt.0.0) then
+ if (wqq_scpho(itypi).ne.0.0) then
Qij=wqq_scpho(itypi)/eps_in
+ alpha_sco=1.d0/alphi_scpho(itypi)
! Qij=0.0
- Ecl = (332.0d0 * Qij) / Rhead
+ Ecl = (332.0d0 * Qij*dexp(-Rhead*alpha_sco)) / Rhead
!c! derivative of Ecl is Gcl...
- dGCLdR = (-332.0d0 * Qij ) / Rhead_sq
+ dGCLdR = (-332.0d0 * Qij*dexp(-Rhead*alpha_sco)* &
+ (Rhead*alpha_sco+1) ) / Rhead_sq
+ if (energy_dec) write(iout,*) "ECL",ECL,Rhead,1.0/rij
else if (wqdip_scpho(2,itypi).gt.0.0d0) then
w1 = wqdip_scpho(1,itypi)
w2 = wqdip_scpho(2,itypi)
Ecl = sparrow / Rhead**2.0d0 &
- hawk / Rhead**4.0d0
!c!-------------------------------------------------------------------
+ if (energy_dec) write(iout,*) "ECLdipdip",ECL,Rhead,&
+ 1.0/rij,sparrow
+
!c! derivative of ecl is Gcl
!c! dF/dr part
dGCLdR = - 2.0d0 * sparrow / Rhead**3.0d0 &
END DO
! if (i.eq.3) print *,i,j,evdwij,epol,Fcav,ECL
+ if (energy_dec) write (iout,'(a22,2i5,4f8.3,f16.3)'), &
+ "escpho:evdw,pol,cav,CL",i,j,evdwij,epol,Fcav,ECL,escpho
escpho=escpho+evdwij+epol+Fcav+ECL
call sc_grad_scpho
enddo
END DO
RETURN
END SUBROUTINE sc_grad_scpho
+ subroutine eprot_pep_phosphate(epeppho)
+ use calc_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.GEO'
+! include 'COMMON.VAR'
+! include 'COMMON.LOCAL'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.NAMES'
+! include 'COMMON.INTERACT'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CALC'
+! include 'COMMON.CONTROL'
+! include 'COMMON.SBRIDGE'
+ logical :: lprn
+!el local variables
+ integer :: iint,itypi,itypi1,itypj,subchap
+ real(kind=8) :: rrij,xi,yi,zi,sig,rij_shift,fac,e1,e2,sigm,epsi
+ real(kind=8) :: evdw,sig0ij
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init,aa,bb,ssgradlipi,ssgradlipj, &
+ sslipi,sslipj,faclip
+ integer :: ii
+ real(kind=8) :: fracinbuf
+ real (kind=8) :: epeppho
+ real (kind=8),dimension(4):: ener
+ real(kind=8) :: b1,b2,b3,b4,egb,eps_in,eps_inout_fac,eps_out
+ real(kind=8) :: ECL,Elj,Equad,Epol,eheadtail,rhead,dGCLOM2,&
+ sqom1,sqom2,sqom12,c1,c2,pom,Lambf,sparrow,&
+ Chif,ChiLambf,bat,eagle,top,bot,botsq,Fcav,dtop,dFdR,dFdOM1,&
+ dFdOM2,w1,w2,dGCLdR,dFdL,dFdOM12,dbot ,&
+ r1,eps_head,alphapol1,pis,facd2,d2,facd1,d1,erdxj,erdxi,federmaus,&
+ dPOLdR1,dFGBdOM2,dFGBdR1,dPOLdFGB1,RR1,MomoFac1,hawk,d1i,d1j,&
+ sig1,sig2,chis12,chis2,ee1,fgb1,a12sq,chis1,Rhead_sq,Qij,dFGBdOM1
+ real(kind=8),dimension(3,2)::chead,erhead_tail
+ real(kind=8),dimension(3) :: Rhead_distance,ertail,erhead
+ integer troll
+ real (kind=8) :: dcosom1(3),dcosom2(3)
+ epeppho=0.0d0
+! do i=1,nres_molec(1)
+ do i=ibond_start,ibond_end
+ if (itype(i,1).eq.ntyp1_molec(1)) cycle
+ itypi = itype(i,1)
+ dsci_inv = vbld_inv(i+1)/2.0
+ dxi = dc_norm(1,i)
+ dyi = dc_norm(2,i)
+ dzi = dc_norm(3,i)
+ xi=(c(1,i)+c(1,i+1))/2.0
+ yi=(c(2,i)+c(2,i+1))/2.0
+ zi=(c(3,i)+c(3,i+1))/2.0
+ 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=nres_molec(1)+1,nres_molec(2)+nres_molec(1)-1
+ itypj= itype(j,2)
+ if ((itype(j,2).eq.ntyp1_molec(2)).or.&
+ (itype(j+1,2).eq.ntyp1_molec(2))) cycle
+ xj=(c(1,j)+c(1,j+1))/2.0
+ yj=(c(2,j)+c(2,j+1))/2.0
+ zj=(c(3,j)+c(3,j+1))/2.0
+ 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
+ rrij = 1.0D0 / ( xj*xj + yj*yj + zj*zj)
+ rij = dsqrt(rrij)
+ dxj = dc_norm( 1,j )
+ dyj = dc_norm( 2,j )
+ dzj = dc_norm( 3,j )
+ dscj_inv = vbld_inv(j+1)/2.0
+! Gay-berne var's
+ sig0ij = sigma_peppho
+ chi1=0.0d0
+ chi2=0.0d0
+ chi12 = chi1 * chi2
+ chip1=0.0d0
+ chip2=0.0d0
+ chip12 = chip1 * chip2
+ chis1 = 0.0d0
+ chis2 = 0.0d0
+ chis12 = chis1 * chis2
+ sig1 = sigmap1_peppho
+ sig2 = sigmap2_peppho
+! write (*,*) "sig1 = ", sig1
+! write (*,*) "sig1 = ", sig1
+! write (*,*) "sig2 = ", sig2
+! alpha factors from Fcav/Gcav
+ alf1 = 0.0d0
+ alf2 = 0.0d0
+ alf12 = 0.0d0
+ b1 = alphasur_peppho(1)
+! b1=0.0d0
+ b2 = alphasur_peppho(2)
+ b3 = alphasur_peppho(3)
+ b4 = alphasur_peppho(4)
+ CALL sc_angular
+ sqom1=om1*om1
+ evdwij = 0.0d0
+ ECL = 0.0d0
+ Elj = 0.0d0
+ Equad = 0.0d0
+ Epol = 0.0d0
+ Fcav=0.0d0
+ eheadtail = 0.0d0
+ dGCLdR=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
+ rij_shift = rij
+ fac = rij_shift**expon
+ c1 = fac * fac * aa_peppho
+! c1 = 0.0d0
+ c2 = fac * bb_peppho
+! c2 = 0.0d0
+ evdwij = c1 + c2
+! Now cavity....................
+ eagle = dsqrt(1.0/rij_shift)
+ top = b1 * ( eagle + b2 * 1.0/rij_shift - b3 )
+ bot = 1.0d0 + b4 * (1.0/rij_shift ** 12.0d0)
+ botsq = bot * bot
+ Fcav = top / bot
+ dtop = b1 * ((1.0/ (2.0d0 * eagle)) + (b2))
+ dbot = 12.0d0 * b4 * (1.0/rij_shift) ** 11.0d0
+ dFdR = ((dtop * bot - top * dbot) / botsq)
+ w1 = wqdip_peppho(1)
+ w2 = wqdip_peppho(2)
+! w1=0.0d0
+! w2=0.0d0
+! pis = sig0head_scbase(itypi,itypj)
+! eps_head = epshead_scbase(itypi,itypj)
+!c!-------------------------------------------------------------------
+
+!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 * om1
+ hawk = w2 * (1.0d0 - sqom1)
+ Ecl = sparrow * rij_shift**2.0d0 &
+ - hawk * rij_shift**4.0d0
+!c!-------------------------------------------------------------------
+!c! derivative of ecl is Gcl
+!c! dF/dr part
+! rij_shift=5.0
+ dGCLdR = - 2.0d0 * sparrow * rij_shift**3.0d0 &
+ + 4.0d0 * hawk * rij_shift**5.0d0
+!c! dF/dom1
+ dGCLdOM1 = (w1) * (rij_shift**2.0d0)
+!c! dF/dom2
+ dGCLdOM2 = (2.0d0 * w2 * om1) * (rij_shift ** 4.0d0)
+ eom1 = dGCLdOM1+dGCLdOM2
+ eom2 = 0.0
+
+ fac = -expon * (c1 + evdwij) * rij_shift+dFdR+dGCLdR
+! fac=0.0
+ gg(1) = fac*xj*rij
+ gg(2) = fac*yj*rij
+ gg(3) = fac*zj*rij
+ do k=1,3
+ gvdwc_peppho(k,j) = gvdwc_peppho(k,j) +gg(k)/2.0
+ gvdwc_peppho(k,j+1) = gvdwc_peppho(k,j+1) +gg(k)/2.0
+ gvdwc_peppho(k,i) = gvdwc_peppho(k,i) -gg(k)/2.0
+ gvdwc_peppho(k,i+1) = gvdwc_peppho(k,i+1) -gg(k)/2.0
+ gg(k)=0.0
+ enddo
+ DO k = 1, 3
+ dcosom1(k) = rij* (dc_norm(k,i) - om1 * erij(k))
+ dcosom2(k) = rij* (dc_norm(k,j) - om2 * erij(k))
+ gg(k) = gg(k) + eom1 * dcosom1(k)! + eom2 * dcosom2(k)
+ gvdwc_peppho(k,j)= gvdwc_peppho(k,j) +0.5*( gg(k)) !&
+! - (eom2*(erij(k)-om2*dc_norm(k,j)))*dscj_inv*2.0
+ gvdwc_peppho(k,j+1)= gvdwc_peppho(k,j+1) +0.5*( gg(k)) !&
+! + (eom2*(erij(k)-om2*dc_norm(k,j)))*dscj_inv*2.0
+ gvdwc_peppho(k,i)= gvdwc_peppho(k,i) -0.5*( gg(k)) &
+ - (eom1*(erij(k)-om1*dc_norm(k,i)))*dsci_inv*2.0
+ gvdwc_peppho(k,i+1)= gvdwc_peppho(k,i+1) - 0.5*( gg(k)) &
+ + (eom1*(erij(k)-om1*dc_norm(k,i)))*dsci_inv*2.0
+ enddo
+ epeppho=epeppho+evdwij+Fcav+ECL
+! print *,i,j,evdwij,Fcav,ECL,rij_shift
+ enddo
+ enddo
+ end subroutine eprot_pep_phosphate
end module energy
projects / unres4.git / blobdiff