X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD%2Fenergy_p_new_barrier.F;h=7cf67dba1c3afaa0bc6a74f506e8e8fd5c272915;hb=a4e5a256e77beab0f13a18e273b46707d587060a;hp=a2eb0f9288190efdc4d89a714dd53ea4ad1afcad;hpb=29858260fdd1ede4a6a430ee7883f998abaac045;p=unres.git diff --git a/source/unres/src_MD/energy_p_new_barrier.F b/source/unres/src_MD/energy_p_new_barrier.F index a2eb0f9..7cf67db 100644 --- a/source/unres/src_MD/energy_p_new_barrier.F +++ b/source/unres/src_MD/energy_p_new_barrier.F @@ -27,6 +27,7 @@ cMS$ATTRIBUTES C :: proc_proc #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() @@ -131,6 +132,31 @@ C 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 @@ -228,6 +254,15 @@ cd print *,'nterm=',nterm 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 @@ -329,6 +364,11 @@ C 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 @@ -426,20 +466,29 @@ cMS$ATTRIBUTES C :: proc_proc 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 @@ -546,7 +595,11 @@ c enddo & 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 @@ -560,7 +613,11 @@ c enddo & 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 @@ -576,7 +633,11 @@ c enddo & 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 @@ -1048,6 +1109,13 @@ C------------------------------------------------------------------------ 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, @@ -1055,8 +1123,9 @@ C------------------------------------------------------------------------ & 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)'/ @@ -1078,8 +1147,13 @@ C------------------------------------------------------------------------ & '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, @@ -1088,7 +1162,9 @@ C------------------------------------------------------------------------ & 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)'/ @@ -1109,8 +1185,13 @@ C------------------------------------------------------------------------ & '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 @@ -4637,7 +4718,7 @@ C Derivatives of the "mean" values in gamma1 and gamma2. & '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 @@ -4928,13 +5009,15 @@ C 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 @@ -5758,6 +5841,15 @@ C Proline-Proline pair is a special case... 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 @@ -5858,6 +5950,187 @@ cd write (iout,*) 'edihcnstr',edihcnstr 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)