#ifdef MPI
c print*,"ETOTAL Processor",fg_rank," absolute rank",myrank,
c & " nfgtasks",nfgtasks
+ call flush(iout)
if (nfgtasks.gt.1) then
#ifdef MPI
time00=MPI_Wtime()
C Calculate electrostatic (H-bonding) energy of the main chain.
C
107 continue
+C BARTEK for dfa test!
+ if (wdfa_dist.gt.0) then
+ call edfad(edfadis)
+ else
+ edfadis=0
+ endif
+c print*, 'edfad is finished!', edfadis
+ if (wdfa_tor.gt.0) then
+ call edfat(edfator)
+ else
+ edfator=0
+ endif
+c print*, 'edfat is finished!', edfator
+ if (wdfa_nei.gt.0) then
+ call edfan(edfanei)
+ else
+ edfanei=0
+ endif
+c print*, 'edfan is finished!', edfanei
+ if (wdfa_beta.gt.0) then
+ call edfab(edfabet)
+ else
+ edfabet=0
+ endif
+c print*, 'edfab is finished!', edfabet
cmc
cmc Sep-06: egb takes care of dynamic ss bonds too
cmc
- if (dyn_ss) call dyn_set_nss
+c if (dyn_ss) call dyn_set_nss
c print *,"Processor",myrank," computed USCSC"
#ifdef TIMING
etors=0
edihcnstr=0
endif
+
+ if (constr_homology.ge.1) then
+ call e_modeller(ehomology_constr)
+ else
+ ehomology_constr=0.0d0
+ endif
+
+
+c write(iout,*) ehomology_constr
c print *,"Processor",myrank," computed Utor"
C
C 6/23/01 Calculate double-torsional energy
energia(21)=esccor
energia(22)=evdw_p
energia(23)=evdw_m
+ energia(24)=ehomology_constr
+ energia(25)=edfadis
+ energia(26)=edfator
+ energia(27)=edfanei
+ energia(28)=edfabet
c print *," Processor",myrank," calls SUM_ENERGY"
call sum_energy(energia,.true.)
+ if (dyn_ss) call dyn_set_nss
c print *," Processor",myrank," left SUM_ENERGY"
#ifdef TIMING
#ifdef MPI
estr=energia(17)
Uconst=energia(20)
esccor=energia(21)
+ ehomology_constr=energia(24)
+ edfadis=energia(25)
+ edfator=energia(26)
+ edfanei=energia(27)
+ edfabet=energia(28)
#ifdef SPLITELE
etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1
& +wang*ebe+wtor*etors+wscloc*escloc
& +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5
& +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
& +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
- & +wbond*estr+Uconst+wsccor*esccor
+ & +wbond*estr+Uconst+wsccor*esccor+ehomology_constr
+ & +wdfa_dist*edfadis+wdfa_tor*edfator+wdfa_nei*edfanei
+ & +wdfa_beta*edfabet
#else
etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1)
& +wang*ebe+wtor*etors+wscloc*escloc
& +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5
& +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
& +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
- & +wbond*estr+Uconst+wsccor*esccor
+ & +wbond*estr+Uconst+wsccor*esccor+ehomology_constr
+ & +wdfa_dist*edfadis+wdfa_tor*edfator+wdfa_nei*edfanei
+ & +wdfa_beta*edfabet
#endif
energia(0)=etot
c detecting NaNQ
& wcorr5*gradcorr5_long(j,i)+
& wcorr6*gradcorr6_long(j,i)+
& wturn6*gcorr6_turn_long(j,i)+
- & wstrain*ghpbc(j,i)
+ & wstrain*ghpbc(j,i)+
+ & wdfa_dist*gdfad(j,i)+
+ & wdfa_tor*gdfat(j,i)+
+ & wdfa_nei*gdfan(j,i)+
+ & wdfa_beta*gdfab(j,i)
enddo
enddo
#else
& wcorr5*gradcorr5_long(j,i)+
& wcorr6*gradcorr6_long(j,i)+
& wturn6*gcorr6_turn_long(j,i)+
- & wstrain*ghpbc(j,i)
+ & wstrain*ghpbc(j,i)+
+ & wdfa_dist*gdfad(j,i)+
+ & wdfa_tor*gdfat(j,i)+
+ & wdfa_nei*gdfan(j,i)+
+ & wdfa_beta*gdfab(j,i)
enddo
enddo
#endif
& wcorr5*gradcorr5_long(j,i)+
& wcorr6*gradcorr6_long(j,i)+
& wturn6*gcorr6_turn_long(j,i)+
- & wstrain*ghpbc(j,i)
+ & wstrain*ghpbc(j,i)+
+ & wdfa_dist*gdfad(j,i)+
+ & wdfa_tor*gdfat(j,i)+
+ & wdfa_nei*gdfan(j,i)+
+ & wdfa_beta*gdfab(j,i)
enddo
enddo
#endif
estr=energia(17)
Uconst=energia(20)
esccor=energia(21)
+ ehomology_constr=energia(24)
+C Bartek
+ edfadis = energia(25)
+ edfator = energia(26)
+ edfanei = energia(27)
+ edfabet = energia(28)
+
#ifdef SPLITELE
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,
& estr,wbond,ebe,wang,
& ecorr,wcorr,
& ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
& eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,
- & edihcnstr,ebr*nss,
- & Uconst,etot
+ & edihcnstr,ehomology_constr, ebr*nss,
+ & Uconst,edfadis,wdfa_dist,edfator,wdfa_tor,edfanei,wdfa_nei,
+ & edfabet,wdfa_beta,etot
10 format (/'Virtual-chain energies:'//
& 'EVDW= ',1pE16.6,' WEIGHT=',1pE16.6,' (SC-SC)'/
& 'EVDW2= ',1pE16.6,' WEIGHT=',1pE16.6,' (SC-p)'/
& 'ETURN6=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 6th order)'/
& 'ESCCOR=',1pE16.6,' WEIGHT=',1pE16.6,' (backbone-rotamer corr)'/
& 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
+ & 'H_CONS=',1pE16.6,' (Homology model constraints energy)'/
& 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
& 'UCONST= ',1pE16.6,' (Constraint energy)'/
+ & 'EDFAD= ',1pE16.6,' WEIGHT=',1pE16.6,' (DFA distance energy)'/
+ & 'EDFAT= ',1pE16.6,' WEIGHT=',1pE16.6,' (DFA torsion energy)'/
+ & 'EDFAN= ',1pE16.6,' WEIGHT=',1pE16.6,' (DFA NCa energy)'/
+ & 'EDFAB= ',1pE16.6,' WEIGHT=',1pE16.6,' (DFA Beta energy)'/
& 'ETOT= ',1pE16.6,' (total)')
#else
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,
& ecorr,wcorr,
& ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
& eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr,
- & ebr*nss,Uconst,etot
+ & ehomology_constr,ebr*nss,Uconst,edfadis,wdfa_dist,edfator,
+ & wdfa_tor,edfanei,wdfa_nei,edfabet,wdfa_beta,
+ & etot
10 format (/'Virtual-chain energies:'//
& 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
& 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
& 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/
& 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/
& 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
+ & 'H_CONS=',1pE16.6,' (Homology model constraints energy)'/
& 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
& 'UCONST=',1pE16.6,' (Constraint energy)'/
+ & 'EDFAD= ',1pE16.6,' WEIGHT=',1pD16.6,' (DFA distance energy)'/
+ & 'EDFAT= ',1pE16.6,' WEIGHT=',1pD16.6,' (DFA torsion energy)'/
+ & 'EDFAN= ',1pE16.6,' WEIGHT=',1pD16.6,' (DFA NCa energy)'/
+ & 'EDFAB= ',1pE16.6,' WEIGHT=',1pD16.6,' (DFA Beta energy)'/
& 'ETOT= ',1pE16.6,' (total)')
#endif
return
& 'ebend',i,ethetai
if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1
if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2
- gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)
+ gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)+gloc(nphi+i-2,icg)
enddo
C Ufff.... We've done all this!!!
return
enddo
enddo
10 continue
- if (lprn1) write (iout,'(i2,3f8.1,9h ethetai ,f10.5)')
- & i,theta(i)*rad2deg,phii*rad2deg,
+c lprn1=.true.
+ if (lprn1) write (iout,'(a4,i2,3f8.1,9h ethetai ,f10.5)')
+ & 'ebe', i,theta(i)*rad2deg,phii*rad2deg,
& phii1*rad2deg,ethetai
+c lprn1=.false.
etheta=etheta+ethetai
if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii
if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1
- gloc(nphi+i-2,icg)=wang*dethetai
+ gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*dethetai
enddo
return
end
return
end
c------------------------------------------------------------------------------
+c LICZENIE WIEZOW Z ROWNANIA ENERGII MODELLERA
+ subroutine e_modeller(ehomology_constr)
+ ehomology_constr=0.0
+ write (iout,*) "!!!!!UWAGA, JESTEM W DZIWNEJ PETLI, TEST!!!!!"
+ return
+ end
+C !!!!!!!! NIE CZYTANE !!!!!!!!!!!
+
+c------------------------------------------------------------------------------
subroutine etor_d(etors_d)
etors_d=0.0d0
return
return
end
c----------------------------------------------------------------------------
+c MODELLER restraint function
+ subroutine e_modeller(ehomology_constr)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+
+ integer nnn, i, j, k, ki, irec, l
+ integer katy, odleglosci, test7
+ real*8 odleg, odleg2, odleg3, kat, kat2, kat3, gdih(max_template)
+ real*8 distance(max_template),distancek(max_template),
+ & min_odl,godl(max_template),dih_diff(max_template)
+
+ include 'COMMON.SBRIDGE'
+ include 'COMMON.CHAIN'
+ include 'COMMON.GEO'
+ include 'COMMON.DERIV'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.VAR'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.MD'
+ include 'COMMON.CONTROL'
+
+
+ do i=1,19
+ distancek(i)=9999999.9
+ enddo
+
+
+ odleg=0.0d0
+
+c Pseudo-energy and gradient from homology restraints (MODELLER-like
+c function)
+C AL 5/2/14 - Introduce list of restraints
+ do ii = link_start_homo,link_end_homo
+ i = ires_homo(ii)
+ j = jres_homo(ii)
+ dij=dist(i,j)
+ do k=1,constr_homology
+ distance(k)=odl(k,ii)-dij
+ distancek(k)=
+ & 0.5d0*waga_dist(iset)*distance(k)**2*sigma_odl(k,ii)
+ enddo
+
+ min_odl=minval(distancek)
+#ifdef DEBUG
+ write (iout,*) "ij dij",i,j,dij
+ write (iout,*) "distance",(distance(k),k=1,constr_homology)
+ write (iout,*) "distancek",(distancek(k),k=1,constr_homology)
+ write (iout,* )"min_odl",min_odl
+#endif
+ odleg2=0.0d0
+ do k=1,constr_homology
+c Nie wiem po co to liczycie jeszcze raz!
+c odleg3=-waga_dist(iset)*((distance(i,j,k)**2)/
+c & (2*(sigma_odl(i,j,k))**2))
+ godl(k)=dexp(-distancek(k)+min_odl)
+ odleg2=odleg2+godl(k)
+
+ccc write(iout,779) i,j,k, "odleg2=",odleg2, "odleg3=", odleg3,
+ccc & "dEXP(odleg3)=", dEXP(odleg3),"distance(i,j,k)^2=",
+ccc & distance(i,j,k)**2, "dist(i+1,j+1)=", dist(i+1,j+1),
+ccc & "sigma_odl(i,j,k)=", sigma_odl(i,j,k)
+
+ enddo
+#ifdef DEBUG
+ write (iout,*) "godl",(godl(k),k=1,constr_homology)
+ write (iout,*) "ii i j",ii,i,j," odleg2",odleg2
+#endif
+ odleg=odleg-dLOG(odleg2/constr_homology)+min_odl
+c Gradient
+ sum_godl=odleg2
+ sum_sgodl=0.0
+ do k=1,constr_homology
+c godl=dexp(((-(distance(i,j,k)**2)/(2*(sigma_odl(i,j,k))**2))
+c & *waga_dist(iset))+min_odl
+ sgodl=-godl(k)*distance(k)*sigma_odl(k,ii)*waga_dist(iset)
+ sum_sgodl=sum_sgodl+sgodl
+
+c sgodl2=sgodl2+sgodl
+c write(iout,*) i, j, k, distance(i,j,k), "W GRADIENCIE1"
+c write(iout,*) "constr_homology=",constr_homology
+c write(iout,*) i, j, k, "TEST K"
+ enddo
+
+ grad_odl3=sum_sgodl/(sum_godl*dij)
+
+
+c write(iout,*) i, j, k, distance(i,j,k), "W GRADIENCIE2"
+c write(iout,*) (distance(i,j,k)**2), (2*(sigma_odl(i,j,k))**2),
+c & (-(distance(i,j,k)**2)/(2*(sigma_odl(i,j,k))**2))
+
+ccc write(iout,*) godl, sgodl, grad_odl3
+
+c grad_odl=grad_odl+grad_odl3
+
+ do jik=1,3
+ ggodl=grad_odl3*(c(jik,i)-c(jik,j))
+ccc write(iout,*) c(jik,i+1), c(jik,j+1), (c(jik,i+1)-c(jik,j+1))
+ccc write(iout,746) "GRAD_ODL_1", i, j, jik, ggodl,
+ccc & ghpbc(jik,i+1), ghpbc(jik,j+1)
+ ghpbc(jik,i)=ghpbc(jik,i)+ggodl
+ ghpbc(jik,j)=ghpbc(jik,j)-ggodl
+ccc write(iout,746) "GRAD_ODL_2", i, j, jik, ggodl,
+ccc & ghpbc(jik,i+1), ghpbc(jik,j+1)
+
+ enddo
+ccc write(iout,778)"TEST: odleg2=", odleg2, "DLOG(odleg2)=",
+ccc & dLOG(odleg2),"-odleg=", -odleg
+
+ enddo ! ii
+c Pseudo-energy and gradient from dihedral-angle restraints from
+c homology templates
+c write (iout,*) "End of distance loop"
+c call flush(iout)
+ kat=0.0d0
+c write (iout,*) idihconstr_start_homo,idihconstr_end_homo
+ do i=idihconstr_start_homo,idihconstr_end_homo
+ kat2=0.0d0
+c betai=beta(i,i+1,i+2,i+3)
+ betai = phi(i+3)
+ do k=1,constr_homology
+ dih_diff(k)=pinorm(dih(k,i)-betai)
+c if (dih_diff(i,k).gt.3.14159) dih_diff(i,k)=
+c & -(6.28318-dih_diff(i,k))
+c if (dih_diff(i,k).lt.-3.14159) dih_diff(i,k)=
+c & 6.28318+dih_diff(i,k)
+
+ kat3=-0.5d0*waga_angle(iset)*dih_diff(k)**2*sigma_dih(k,i)
+ gdih(k)=dexp(kat3)
+ kat2=kat2+gdih(k)
+c write(iout,*) "kat2=", kat2, "exp(kat3)=", exp(kat3)
+c write(*,*)""
+ enddo
+#ifdef DEBUG
+ write (iout,*) "i",i," betai",betai," kat2",kat2
+ write (iout,*) "gdih",(gdih(k),k=1,constr_homology)
+#endif
+ if (kat2.le.1.0d-14) cycle
+ kat=kat-dLOG(kat2/constr_homology)
+
+ccc write(iout,778)"TEST: kat2=", kat2, "DLOG(kat2)=",
+ccc & dLOG(kat2), "-kat=", -kat
+
+c ----------------------------------------------------------------------
+c Gradient
+c ----------------------------------------------------------------------
+
+ sum_gdih=kat2
+ sum_sgdih=0.0
+ do k=1,constr_homology
+ sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i)*waga_angle(iset)
+ sum_sgdih=sum_sgdih+sgdih
+ enddo
+ grad_dih3=sum_sgdih/sum_gdih
+
+c write(iout,*)i,k,gdih,sgdih,beta(i+1,i+2,i+3,i+4),grad_dih3
+ccc write(iout,747) "GRAD_KAT_1", i, nphi, icg, grad_dih3,
+ccc & gloc(nphi+i-3,icg)
+ gloc(i,icg)=gloc(i,icg)+grad_dih3
+ccc write(iout,747) "GRAD_KAT_2", i, nphi, icg, grad_dih3,
+ccc & gloc(nphi+i-3,icg)
+
+ enddo
+
+
+c Total energy from homology restraints
+#ifdef DEBUG
+ write (iout,*) "odleg",odleg," kat",kat
+#endif
+ ehomology_constr=odleg+kat
+ return
+
+ 748 format(a8,f12.3,a6,f12.3,a7,f12.3)
+ 747 format(a12,i4,i4,i4,f8.3,f8.3)
+ 746 format(a12,i4,i4,i4,f8.3,f8.3,f8.3)
+ 778 format(a7,1X,f10.3,1X,a4,1X,f10.3,1X,a5,1X,f10.3)
+ 779 format(i3,1X,i3,1X,i2,1X,a7,1X,f7.3,1X,a7,1X,f7.3,1X,a13,1X,
+ & f7.3,1X,a17,1X,f9.3,1X,a10,1X,f8.3,1X,a10,1X,f8.3)
+ end
+
+c------------------------------------------------------------------------------
subroutine etor_d(etors_d)
C 6/23/01 Compute double torsional energy
implicit real*8 (a-h,o-z)