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
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
difi=phii-phi0(i)
if (difi.gt.drange(i)) then
difi=difi-drange(i)
- edihcnstr=edihcnstr+0.25d0*ftors*difi**4
- gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
+ edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
+ gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
else if (difi.lt.-drange(i)) then
difi=difi+drange(i)
- edihcnstr=edihcnstr+0.25d0*ftors*difi**4
- gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
+ edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
+ gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
endif
! write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii,
! & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg)
difi=pinorm(phii-phi0(i))
if (difi.gt.drange(i)) then
difi=difi-drange(i)
- edihcnstr=edihcnstr+0.25d0*ftors*difi**4
- gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
+ edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
+ gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
else if (difi.lt.-drange(i)) then
difi=difi+drange(i)
- edihcnstr=edihcnstr+0.25d0*ftors*difi**4
- gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
+ edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
+ gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
else
difi=0.0
endif
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(3,itypi,itypj)*2.0
+ 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(3,itypj)*2.0
+ 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
!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
integer troll
real (kind=8) :: dcosom1(3),dcosom2(3)
epeppho=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)
dsci_inv = vbld_inv(i+1)/2.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
projects / unres4.git / blobdiff