X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fenergy.F90;h=ddc9833cd43ff65c32bf3957f2752600f8acda9e;hb=bc23440fbe68672d430f71f22f46b11265f003db;hp=677e9705eca6725ba0f2a7b6c50ffc3419ee2484;hpb=58cb5b8fe97f65aba46ef7c9a0d6ab38061d8cab;p=unres4.git diff --git a/source/unres/energy.F90 b/source/unres/energy.F90 index 677e970..ddc9833 100644 --- a/source/unres/energy.F90 +++ b/source/unres/energy.F90 @@ -241,7 +241,7 @@ real(kind=8) :: evdw,evdw1,evdw2,evdw2_14,escloc,ees,eel_loc real(kind=8) :: eello_turn3,eello_turn4,estr,ebe,eliptran,etube, & Eafmforce,ethetacnstr - real(kind=8) :: ecorr,ecorr5,ecorr6,eturn6 + real(kind=8) :: ecorr,ecorr5,ecorr6,eturn6,ehomology_constr ! now energies for nulceic alone parameters real(kind=8) :: evdwpp,eespp,evdwpsb,eelpsb,evdwsb,eelsb,estr_nucl,& ebe_nucl,esbloc,etors_nucl,etors_d_nucl,ecorr_nucl,& @@ -285,7 +285,7 @@ ! allocate(ishield_listbuf(nres)) ! allocate(shield_listbuf(maxcontsshi,nres)) ! endif - +! print *,"wstrain check", wstrain ! print*,"ETOTAL Processor",fg_rank," absolute rank",myrank, ! & " nfgtasks",nfgtasks if (nfgtasks.gt.1) then @@ -745,7 +745,16 @@ !c print *,"Processor",myrank," computed Utor" ! print *,"Processor",myrank," computed Utor" - + if (constr_homology.ge.1) then + call e_modeller(ehomology_constr) +! print *,'iset=',iset,'me=',me,ehomology_constr, +! & 'Processor',fg_rank,' CG group',kolor, +! & ' absolute rank',MyRank +! print *,"tu" + else + ehomology_constr=0.0d0 + endif + ! ! 6/23/01 Calculate double-torsional energy ! @@ -796,8 +805,8 @@ ! ! If performing constraint dynamics, call the constraint energy ! after the equilibration time - if(usampl.and.totT.gt.eq_time) then -!elwrite(iout,*) "afeter multibody hb" + if((usampl).and.(totT.gt.eq_time)) then + write(iout,*) "usampl",usampl call EconstrQ !elwrite(iout,*) "afeter multibody hb" call Econstr_back @@ -966,6 +975,7 @@ energia(49)=epeppho ! energia(50)=ecations_prot_amber energia(50)=ecation_nucl + energia(51)=ehomology_constr call sum_energy(energia,.true.) if (dyn_ss) call dyn_set_nss ! print *," Processor",myrank," left SUM_ENERGY" @@ -1007,7 +1017,7 @@ eliptran,etube, Eafmforce,ethetacnstr real(kind=8) :: evdwpp,eespp,evdwpsb,eelpsb,evdwsb,eelsb,estr_nucl,& ebe_nucl,esbloc,etors_nucl,etors_d_nucl,ecorr_nucl,& - ecorr3_nucl + ecorr3_nucl,ehomology_constr real(kind=8) :: ecation_prot,ecationcation,ecations_prot_amber,& ecation_nucl real(kind=8) :: escbase,epepbase,escpho,epeppho @@ -1094,6 +1104,7 @@ escpho=energia(48) epeppho=energia(49) ecation_nucl=energia(50) + ehomology_constr=energia(51) ! ecations_prot_amber=energia(50) ! energia(41)=ecation_prot @@ -1107,7 +1118,7 @@ +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+wtube*etube& - +Eafmforce+ethetacnstr & + +Eafmforce+ethetacnstr+ehomology_constr & +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl& +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb& +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl& @@ -1121,7 +1132,7 @@ +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+wtube*etube& - +Eafmforce+ethetacnstr & + +Eafmforce+ethetacnstr+ehomology_constr & +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl& +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb& +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl& @@ -1259,7 +1270,7 @@ etube,ethetacnstr,Eafmforce real(kind=8) :: evdwpp,eespp,evdwpsb,eelpsb,evdwsb,eelsb,estr_nucl,& ebe_nucl,esbloc,etors_nucl,etors_d_nucl,ecorr_nucl,& - ecorr3_nucl + ecorr3_nucl,ehomology_constr real(kind=8) :: ecation_prot,ecationcation,ecations_prot_amber,& ecation_nucl real(kind=8) :: escbase,epepbase,escpho,epeppho @@ -1316,6 +1327,8 @@ escpho=energia(48) epeppho=energia(49) ecation_nucl=energia(50) + ehomology_constr=energia(51) + ! ecations_prot_amber=energia(50) #ifdef SPLITELE write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,& @@ -1332,7 +1345,7 @@ 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,& - ecation_nucl,wcatnucl,etot + ecation_nucl,wcatnucl,ehomology_constr,etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ & @@ -1380,6 +1393,7 @@ 'ESCPHO=',1pE16.6,' WEIGHT=',1pD16.6,'(sc-prot nucl-phosphate)'/& 'EPEPPHO=',1pE16.6,' WEIGHT=',1pD16.6,'(pep-prot nucl-phosphate)'/& 'ECATBASE=',1pE16.6,' WEIGHT=',1pD16.6,'(cation nucl-base)'/& + 'H_CONS=',1pE16.6,' (Homology model constraints energy)'/& 'ETOT= ',1pE16.6,' (total)') #else write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,& @@ -1389,14 +1403,14 @@ 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,Eafmforce, & - etube,wtube, & + etube,wtube, ehomology_constr,& 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,& 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,& - ecation_nucl,wcatnucl,etot + ecation_nucl,wcatnucl,ehomology_constr,etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ & @@ -1443,6 +1457,7 @@ 'ESCPHO=',1pE16.6,' WEIGHT=',1pD16.6,'(sc-prot nucl-phosphate)'/& 'EPEPPHO=',1pE16.6,' WEIGHT=',1pD16.6,'(pep-prot nucl-phosphate)'/& 'ECATBASE=',1pE16.6,' WEIGHT=',1pD16.6,'(cation nucl-base)'/& + 'H_CONS=',1pE16.6,' (Homology model constraints energy)'/& 'ETOT= ',1pE16.6,' (total)') #endif return @@ -7205,6 +7220,15 @@ etors_d=0.0d0 return end subroutine etor_d +!----------------------------------------------------------------------------- +c LICZENIE WIEZOW Z ROWNANIA ENERGII MODELLERA + subroutine e_modeller(ehomology_constr) + real(kind=8) :: ehomology_constr + ehomology_constr=0.0d0 + write (iout,*) "!!!!!UWAGA, JESTEM W DZIWNEJ PETLI, TEST!!!!!" + return + end subroutine e_modeller +C !!!!!!!! NIE CZYTANE !!!!!!!!!!! #else !----------------------------------------------------------------------------- subroutine etor(etors) @@ -7551,7 +7575,651 @@ return end subroutine etor_d #endif +!---------------------------------------------------------------------------- +!---------------------------------------------------------------------------- + subroutine e_modeller(ehomology_constr) +! implicit none +! include 'DIMENSIONS' + use MD_data, only: iset + real(kind=8) :: ehomology_constr + integer nnn,i,ii,j,k,ijk,jik,ki,kk,nexl,irec,l + integer katy, odleglosci, test7 + real(kind=8) :: odleg, odleg2, odleg3, kat, kat2, kat3 + real(kind=8) :: Eval,Erot,min_odl + real(kind=8),dimension(constr_homology) :: distance,distancek,godl,dih_diff,gdih, & + gtheta,dscdiff, & + uscdiffk,guscdiff2,guscdiff3,& + theta_diff + + +! +! FP - 30/10/2014 Temporary specifications for homology restraints +! + real(kind=8) :: utheta_i,gutheta_i,sum_gtheta,sum_sgtheta,& + sgtheta + real(kind=8), dimension (nres) :: guscdiff,usc_diff + real(kind=8) :: sum_godl,sgodl,grad_odl3,ggodl,sum_gdih,& + sum_guscdiff,sum_sgdih,sgdih,grad_dih3,usc_diff_i,dxx,dyy,dzz,& + betai,sum_sgodl,dij,max_template +! real(kind=8) :: dist,pinorm +! +! 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' +! include 'COMMON.HOMOLOGY' +! include 'COMMON.QRESTR' +! +! From subroutine Econstr_back +! +! include 'COMMON.NAMES' +! include 'COMMON.TIME1' +! + + + do i=1,max_template + distancek(i)=9999999.9 + enddo + + odleg=0.0d0 + +! Pseudo-energy and gradient from homology restraints (MODELLER-like +! function) +! AL 5/2/14 - Introduce list of restraints +! write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d +#ifdef DEBUG + write(iout,*) "------- dist restrs start -------" +#endif + do ii = link_start_homo,link_end_homo + i = ires_homo(ii) + j = jres_homo(ii) + dij=dist(i,j) +! write (iout,*) "dij(",i,j,") =",dij + nexl=0 + do k=1,constr_homology +! write(iout,*) ii,k,i,j,l_homo(k,ii),dij,odl(k,ii) + if(.not.l_homo(k,ii)) then + nexl=nexl+1 + cycle + endif + distance(k)=odl(k,ii)-dij +! write (iout,*) "distance(",k,") =",distance(k) +! +! For Gaussian-type Urestr +! + distancek(k)=0.5d0*distance(k)**2*sigma_odl(k,ii) ! waga_dist rmvd from Gaussian argument +! write (iout,*) "sigma_odl(",k,ii,") =",sigma_odl(k,ii) +! write (iout,*) "distancek(",k,") =",distancek(k) +! distancek(k)=0.5d0*waga_dist*distance(k)**2*sigma_odl(k,ii) +! +! For Lorentzian-type Urestr +! + if (waga_dist.lt.0.0d0) then + sigma_odlir(k,ii)=dsqrt(1/sigma_odl(k,ii)) + distancek(k)=distance(k)**2/(sigma_odlir(k,ii)* & + (distance(k)**2+sigma_odlir(k,ii)**2)) + endif + enddo + +! min_odl=minval(distancek) + if (nexl.gt.0) then + min_odl=0.0d0 + else + do kk=1,constr_homology + if(l_homo(kk,ii)) then + min_odl=distancek(kk) + exit + endif + enddo + do kk=1,constr_homology + if (l_homo(kk,ii) .and. distancek(kk).lt.min_odl) & + min_odl=distancek(kk) + enddo + endif + +! write (iout,* )"min_odl",min_odl +#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 +#ifdef OLDRESTR + odleg2=0.0d0 +#else + if (waga_dist.ge.0.0d0) then + odleg2=nexl + else + odleg2=0.0d0 + endif +#endif + do k=1,constr_homology +! Nie wiem po co to liczycie jeszcze raz! +! odleg3=-waga_dist(iset)*((distance(i,j,k)**2)/ +! & (2*(sigma_odl(i,j,k))**2)) + if(.not.l_homo(k,ii)) cycle + if (waga_dist.ge.0.0d0) then +! +! For Gaussian-type Urestr +! + godl(k)=dexp(-distancek(k)+min_odl) + odleg2=odleg2+godl(k) +! +! For Lorentzian-type Urestr +! + else + odleg2=odleg2+distancek(k) + endif + +!cc write(iout,779) i,j,k, "odleg2=",odleg2, "odleg3=", odleg3, +!cc & "dEXP(odleg3)=", dEXP(odleg3),"distance(i,j,k)^2=", +!cc & distance(i,j,k)**2, "dist(i+1,j+1)=", dist(i+1,j+1), +!cc & "sigma_odl(i,j,k)=", sigma_odl(i,j,k) + + enddo +! write (iout,*) "godl",(godl(k),k=1,constr_homology) ! exponents +! write (iout,*) "ii i j",ii,i,j," odleg2",odleg2 ! sum of exps +#ifdef DEBUG + write (iout,*) "godl",(godl(k),k=1,constr_homology) ! exponents + write (iout,*) "ii i j",ii,i,j," odleg2",odleg2 ! sum of exps +#endif + if (waga_dist.ge.0.0d0) then +! +! For Gaussian-type Urestr +! + odleg=odleg-dLOG(odleg2/constr_homology)+min_odl +! +! For Lorentzian-type Urestr +! + else + odleg=odleg+odleg2/constr_homology + endif +! +! write (iout,*) "odleg",odleg ! sum of -ln-s +! Gradient +! +! For Gaussian-type Urestr +! + if (waga_dist.ge.0.0d0) sum_godl=odleg2 + sum_sgodl=0.0d0 + do k=1,constr_homology +! godl=dexp(((-(distance(i,j,k)**2)/(2*(sigma_odl(i,j,k))**2)) +! & *waga_dist)+min_odl +! sgodl=-godl(k)*distance(k)*sigma_odl(k,ii)*waga_dist +! + if(.not.l_homo(k,ii)) cycle + if (waga_dist.ge.0.0d0) then +! For Gaussian-type Urestr +! + sgodl=-godl(k)*distance(k)*sigma_odl(k,ii) ! waga_dist rmvd +! +! For Lorentzian-type Urestr +! + else + sgodl=-2*sigma_odlir(k,ii)*(distance(k)/(distance(k)**2+ & + sigma_odlir(k,ii)**2)**2) + endif + sum_sgodl=sum_sgodl+sgodl + +! sgodl2=sgodl2+sgodl +! write(iout,*) i, j, k, distance(i,j,k), "W GRADIENCIE1" +! write(iout,*) "constr_homology=",constr_homology +! write(iout,*) i, j, k, "TEST K" + enddo +! print *, "ok",iset + if (waga_dist.ge.0.0d0) then +! +! For Gaussian-type Urestr +! + grad_odl3=waga_homology(iset)*waga_dist & + *sum_sgodl/(sum_godl*dij) +! print *, "ok" +! +! For Lorentzian-type Urestr +! + else +! Original grad expr modified by analogy w Gaussian-type Urestr grad +! grad_odl3=-waga_homology(iset)*waga_dist*sum_sgodl + grad_odl3=-waga_homology(iset)*waga_dist* & + sum_sgodl/(constr_homology*dij) +! print *, "ok2" + endif +! +! grad_odl3=sum_sgodl/(sum_godl*dij) + + +! write(iout,*) i, j, k, distance(i,j,k), "W GRADIENCIE2" +! write(iout,*) (distance(i,j,k)**2), (2*(sigma_odl(i,j,k))**2), +! & (-(distance(i,j,k)**2)/(2*(sigma_odl(i,j,k))**2)) + +!cc write(iout,*) godl, sgodl, grad_odl3 + +! grad_odl=grad_odl+grad_odl3 + + do jik=1,3 + ggodl=grad_odl3*(c(jik,i)-c(jik,j)) +!cc write(iout,*) c(jik,i+1), c(jik,j+1), (c(jik,i+1)-c(jik,j+1)) +!cc write(iout,746) "GRAD_ODL_1", i, j, jik, ggodl, +!cc & ghpbc(jik,i+1), ghpbc(jik,j+1) + ghpbc(jik,i)=ghpbc(jik,i)+ggodl + ghpbc(jik,j)=ghpbc(jik,j)-ggodl +!cc write(iout,746) "GRAD_ODL_2", i, j, jik, ggodl, +!cc & ghpbc(jik,i+1), ghpbc(jik,j+1) +! if (i.eq.25.and.j.eq.27) then +! write(iout,*) "jik",jik,"i",i,"j",j +! write(iout,*) "sum_sgodl",sum_sgodl,"sgodl",sgodl +! write(iout,*) "grad_odl3",grad_odl3 +! write(iout,*) "c(",jik,i,")",c(jik,i),"c(",jik,j,")",c(jik,j) +! write(iout,*) "ggodl",ggodl +! write(iout,*) "ghpbc(",jik,i,")", +! & ghpbc(jik,i),"ghpbc(",jik,j,")", +! & ghpbc(jik,j) +! endif + enddo +!cc write(iout,778)"TEST: odleg2=", odleg2, "DLOG(odleg2)=", +!cc & dLOG(odleg2),"-odleg=", -odleg + + enddo ! ii-loop for dist +#ifdef DEBUG + write(iout,*) "------- dist restrs end -------" +! if (waga_angle.eq.1.0d0 .or. waga_theta.eq.1.0d0 .or. +! & waga_d.eq.1.0d0) call sum_gradient +#endif +! Pseudo-energy and gradient from dihedral-angle restraints from +! homology templates +! write (iout,*) "End of distance loop" +! call flush(iout) + kat=0.0d0 +! write (iout,*) idihconstr_start_homo,idihconstr_end_homo +#ifdef DEBUG + write(iout,*) "------- dih restrs start -------" + do i=idihconstr_start_homo,idihconstr_end_homo + write (iout,*) "gloc_init(",i,icg,")",gloc(i,icg) + enddo +#endif + do i=idihconstr_start_homo,idihconstr_end_homo + kat2=0.0d0 +! betai=beta(i,i+1,i+2,i+3) + betai = phi(i) +! write (iout,*) "betai =",betai + do k=1,constr_homology + dih_diff(k)=pinorm(dih(k,i)-betai) +!d write (iout,'(a8,2i4,2f15.8)') "dih_diff",i,k,dih_diff(k) +!d & ,sigma_dih(k,i) +! if (dih_diff(i,k).gt.3.14159) dih_diff(i,k)= +! & -(6.28318-dih_diff(i,k)) +! if (dih_diff(i,k).lt.-3.14159) dih_diff(i,k)= +! & 6.28318+dih_diff(i,k) +#ifdef OLD_DIHED + kat3=-0.5d0*dih_diff(k)**2*sigma_dih(k,i) ! waga_angle rmvd from Gaussian argument +#else + kat3=(dcos(dih_diff(k))-1)*sigma_dih(k,i) ! waga_angle rmvd from Gaussian argument +#endif +! kat3=-0.5d0*waga_angle*dih_diff(k)**2*sigma_dih(k,i) + gdih(k)=dexp(kat3) + kat2=kat2+gdih(k) +! write(iout,*) "kat2=", kat2, "exp(kat3)=", exp(kat3) +! write(*,*)"" + enddo +! write (iout,*) "gdih",(gdih(k),k=1,constr_homology) ! exps +! write (iout,*) "i",i," betai",betai," kat2",kat2 ! sum of exps +#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) +! write (iout,*) "kat",kat ! sum of -ln-s + +!cc write(iout,778)"TEST: kat2=", kat2, "DLOG(kat2)=", +!cc & dLOG(kat2), "-kat=", -kat + +! ---------------------------------------------------------------------- +! Gradient +! ---------------------------------------------------------------------- + + sum_gdih=kat2 + sum_sgdih=0.0d0 + do k=1,constr_homology +#ifdef OLD_DIHED + sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i) ! waga_angle rmvd +#else + sgdih=-gdih(k)*dsin(dih_diff(k))*sigma_dih(k,i) ! waga_angle rmvd +#endif +! sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i)*waga_angle + sum_sgdih=sum_sgdih+sgdih + enddo +! grad_dih3=sum_sgdih/sum_gdih + grad_dih3=waga_homology(iset)*waga_angle*sum_sgdih/sum_gdih +! print *, "ok3" + +! write(iout,*)i,k,gdih,sgdih,beta(i+1,i+2,i+3,i+4),grad_dih3 +!cc write(iout,747) "GRAD_KAT_1", i, nphi, icg, grad_dih3, +!cc & gloc(nphi+i-3,icg) + gloc(i-3,icg)=gloc(i-3,icg)+grad_dih3 +! if (i.eq.25) then +! write(iout,*) "i",i,"icg",icg,"gloc(",i,icg,")",gloc(i,icg) +! endif +!cc write(iout,747) "GRAD_KAT_2", i, nphi, icg, grad_dih3, +!cc & gloc(nphi+i-3,icg) + + enddo ! i-loop for dih +#ifdef DEBUG + write(iout,*) "------- dih restrs end -------" +#endif + +! Pseudo-energy and gradient for theta angle restraints from +! homology templates +! FP 01/15 - inserted from econstr_local_test.F, loop structure +! adapted + +! +! For constr_homology reference structures (FP) +! +! Uconst_back_tot=0.0d0 + Eval=0.0d0 + Erot=0.0d0 +! Econstr_back legacy + do i=1,nres +! do i=ithet_start,ithet_end + dutheta(i)=0.0d0 + enddo +! do i=loc_start,loc_end + do i=-1,nres + do j=1,3 + duscdiff(j,i)=0.0d0 + duscdiffx(j,i)=0.0d0 + enddo + enddo +! +! do iref=1,nref +! write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end +! write (iout,*) "waga_theta",waga_theta + if (waga_theta.gt.0.0d0) then +#ifdef DEBUG + write (iout,*) "usampl",usampl + write(iout,*) "------- theta restrs start -------" +! do i=ithet_start,ithet_end +! write (iout,*) "gloc_init(",nphi+i,icg,")",gloc(nphi+i,icg) +! enddo +#endif +! write (iout,*) "maxres",maxres,"nres",nres + + do i=ithet_start,ithet_end +! +! do i=1,nfrag_back +! ii = ifrag_back(2,i,iset)-ifrag_back(1,i,iset) +! +! Deviation of theta angles wrt constr_homology ref structures +! + utheta_i=0.0d0 ! argument of Gaussian for single k + gutheta_i=0.0d0 ! Sum of Gaussians over constr_homology ref structures +! do j=ifrag_back(1,i,iset)+2,ifrag_back(2,i,iset) ! original loop +! over residues in a fragment +! write (iout,*) "theta(",i,")=",theta(i) + do k=1,constr_homology +! +! dtheta_i=theta(j)-thetaref(j,iref) +! dtheta_i=thetaref(k,i)-theta(i) ! original form without indexing + theta_diff(k)=thetatpl(k,i)-theta(i) +!d write (iout,'(a8,2i4,2f15.8)') "theta_diff",i,k,theta_diff(k) +!d & ,sigma_theta(k,i) + +! + utheta_i=-0.5d0*theta_diff(k)**2*sigma_theta(k,i) ! waga_theta rmvd from Gaussian argument +! utheta_i=-0.5d0*waga_theta*theta_diff(k)**2*sigma_theta(k,i) ! waga_theta? + gtheta(k)=dexp(utheta_i) ! + min_utheta_i? + gutheta_i=gutheta_i+gtheta(k) ! Sum of Gaussians (pk) +! Gradient for single Gaussian restraint in subr Econstr_back +! dutheta(j-2)=dutheta(j-2)+wfrag_back(1,i,iset)*dtheta_i/(ii-1) +! + enddo +! write (iout,*) "gtheta",(gtheta(k),k=1,constr_homology) ! exps +! write (iout,*) "i",i," gutheta_i",gutheta_i ! sum of exps + +! +! Gradient for multiple Gaussian restraint + sum_gtheta=gutheta_i + sum_sgtheta=0.0d0 + do k=1,constr_homology +! New generalized expr for multiple Gaussian from Econstr_back + sgtheta=-gtheta(k)*theta_diff(k)*sigma_theta(k,i) ! waga_theta rmvd +! +! sgtheta=-gtheta(k)*theta_diff(k)*sigma_theta(k,i)*waga_theta ! right functional form? + sum_sgtheta=sum_sgtheta+sgtheta ! cum variable + enddo +! Final value of gradient using same var as in Econstr_back + gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg) & + +sum_sgtheta/sum_gtheta*waga_theta & + *waga_homology(iset) +! print *, "ok4" + +! dutheta(i-2)=sum_sgtheta/sum_gtheta*waga_theta +! & *waga_homology(iset) +! dutheta(i)=sum_sgtheta/sum_gtheta +! +! Uconst_back=Uconst_back+waga_theta*utheta(i) ! waga_theta added as weight + Eval=Eval-dLOG(gutheta_i/constr_homology) +! write (iout,*) "utheta(",i,")=",utheta(i) ! -ln of sum of exps +! write (iout,*) "Uconst_back",Uconst_back ! sum of -ln-s +! Uconst_back=Uconst_back+utheta(i) + enddo ! (i-loop for theta) +#ifdef DEBUG + write(iout,*) "------- theta restrs end -------" +#endif + endif +! +! Deviation of local SC geometry +! +! Separation of two i-loops (instructed by AL - 11/3/2014) +! +! write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end +! write (iout,*) "waga_d",waga_d + +#ifdef DEBUG + write(iout,*) "------- SC restrs start -------" + write (iout,*) "Initial duscdiff,duscdiffx" + do i=loc_start,loc_end + write (iout,*) i,(duscdiff(jik,i),jik=1,3), & + (duscdiffx(jik,i),jik=1,3) + enddo +#endif + do i=loc_start,loc_end + usc_diff_i=0.0d0 ! argument of Gaussian for single k + guscdiff(i)=0.0d0 ! Sum of Gaussians over constr_homology ref structures +! do j=ifrag_back(1,i,iset)+1,ifrag_back(2,i,iset)-1 ! Econstr_back legacy +! write(iout,*) "xxtab, yytab, zztab" +! write(iout,'(i5,3f8.2)') i,xxtab(i),yytab(i),zztab(i) + do k=1,constr_homology +! + dxx=-xxtpl(k,i)+xxtab(i) ! Diff b/w x component of ith SC vector in model and kth ref str? +! Original sign inverted for calc of gradients (s. Econstr_back) + dyy=-yytpl(k,i)+yytab(i) ! ibid y + dzz=-zztpl(k,i)+zztab(i) ! ibid z +! write(iout,*) "dxx, dyy, dzz" +!d write(iout,'(2i5,4f8.2)') k,i,dxx,dyy,dzz,sigma_d(k,i) +! + usc_diff_i=-0.5d0*(dxx**2+dyy**2+dzz**2)*sigma_d(k,i) ! waga_d rmvd from Gaussian argument +! usc_diff(i)=-0.5d0*waga_d*(dxx**2+dyy**2+dzz**2)*sigma_d(k,i) ! waga_d? +! uscdiffk(k)=usc_diff(i) + guscdiff2(k)=dexp(usc_diff_i) ! without min_scdiff +! write(iout,*) "i",i," k",k," sigma_d",sigma_d(k,i), +! & " guscdiff2",guscdiff2(k) + guscdiff(i)=guscdiff(i)+guscdiff2(k) !Sum of Gaussians (pk) +! write (iout,'(i5,6f10.5)') j,xxtab(j),yytab(j),zztab(j), +! & xxref(j),yyref(j),zzref(j) + enddo +! +! Gradient +! +! Generalized expression for multiple Gaussian acc to that for a single +! Gaussian in Econstr_back as instructed by AL (FP - 03/11/2014) +! +! Original implementation +! sum_guscdiff=guscdiff(i) +! +! sum_sguscdiff=0.0d0 +! do k=1,constr_homology +! sguscdiff=-guscdiff2(k)*dscdiff(k)*sigma_d(k,i)*waga_d !waga_d? +! sguscdiff=-guscdiff3(k)*dscdiff(k)*sigma_d(k,i)*waga_d ! w min_uscdiff +! sum_sguscdiff=sum_sguscdiff+sguscdiff +! enddo +! +! Implementation of new expressions for gradient (Jan. 2015) +! +! grad_uscdiff=sum_sguscdiff/(sum_guscdiff*dtab) !? + do k=1,constr_homology +! +! New calculation of dxx, dyy, and dzz corrected by AL (07/11), was missing and wrong +! before. Now the drivatives should be correct +! + dxx=-xxtpl(k,i)+xxtab(i) ! Diff b/w x component of ith SC vector in model and kth ref str? +! Original sign inverted for calc of gradients (s. Econstr_back) + dyy=-yytpl(k,i)+yytab(i) ! ibid y + dzz=-zztpl(k,i)+zztab(i) ! ibid z + sum_guscdiff=guscdiff2(k)* &!(dsqrt(dxx*dxx+dyy*dyy+dzz*dzz))* -> wrong! + sigma_d(k,i) ! for the grad wrt r' +! sum_sguscdiff=sum_sguscdiff+sum_guscdiff + +! +! New implementation + sum_guscdiff = waga_homology(iset)*waga_d*sum_guscdiff + do jik=1,3 + duscdiff(jik,i-1)=duscdiff(jik,i-1)+ & + sum_guscdiff*(dXX_C1tab(jik,i)*dxx+ & + dYY_C1tab(jik,i)*dyy+dZZ_C1tab(jik,i)*dzz)/guscdiff(i) + duscdiff(jik,i)=duscdiff(jik,i)+ & + sum_guscdiff*(dXX_Ctab(jik,i)*dxx+ & + dYY_Ctab(jik,i)*dyy+dZZ_Ctab(jik,i)*dzz)/guscdiff(i) + duscdiffx(jik,i)=duscdiffx(jik,i)+ & + sum_guscdiff*(dXX_XYZtab(jik,i)*dxx+ & + dYY_XYZtab(jik,i)*dyy+dZZ_XYZtab(jik,i)*dzz)/guscdiff(i) +! print *, "ok5" +! +#ifdef DEBUG +! write(iout,*) "jik",jik,"i",i + write(iout,*) "dxx, dyy, dzz" + write(iout,'(2i5,3f8.2)') k,i,dxx,dyy,dzz + write(iout,*) "guscdiff2(",k,")",guscdiff2(k) + write(iout,*) "sum_sguscdiff",sum_guscdiff,waga_homology(iset),waga_d + write(iout,*) "dXX_Ctab(",jik,i,")",dXX_Ctab(jik,i) + write(iout,*) "dYY_Ctab(",jik,i,")",dYY_Ctab(jik,i) + write(iout,*) "dZZ_Ctab(",jik,i,")",dZZ_Ctab(jik,i) + write(iout,*) "dXX_C1tab(",jik,i,")",dXX_C1tab(jik,i) + write(iout,*) "dYY_C1tab(",jik,i,")",dYY_C1tab(jik,i) + write(iout,*) "dZZ_C1tab(",jik,i,")",dZZ_C1tab(jik,i) + write(iout,*) "dXX_XYZtab(",jik,i,")",dXX_XYZtab(jik,i) + write(iout,*) "dYY_XYZtab(",jik,i,")",dYY_XYZtab(jik,i) + write(iout,*) "dZZ_XYZtab(",jik,i,")",dZZ_XYZtab(jik,i) + write(iout,*) "duscdiff(",jik,i-1,")",duscdiff(jik,i-1) + write(iout,*) "duscdiff(",jik,i,")",duscdiff(jik,i) + write(iout,*) "duscdiffx(",jik,i,")",duscdiffx(jik,i) +! endif +#endif + enddo + enddo +! print *, "ok6" +! +! uscdiff(i)=-dLOG(guscdiff(i)/(ii-1)) ! Weighting by (ii-1) required? +! usc_diff(i)=-dLOG(guscdiff(i)/constr_homology) ! + min_uscdiff ? +! +! write (iout,*) i," uscdiff",uscdiff(i) +! +! Put together deviations from local geometry + +! Uconst_back=Uconst_back+wfrag_back(1,i,iset)*utheta(i)+ +! & wfrag_back(3,i,iset)*uscdiff(i) + Erot=Erot-dLOG(guscdiff(i)/constr_homology) +! write (iout,*) "usc_diff(",i,")=",usc_diff(i) ! -ln of sum of exps +! write (iout,*) "Uconst_back",Uconst_back ! cum sum of -ln-s +! Uconst_back=Uconst_back+usc_diff(i) +! +! Gradient of multiple Gaussian restraint (FP - 04/11/2014 - right?) +! +! New implment: multiplied by sum_sguscdiff +! + + enddo ! (i-loop for dscdiff) + +! endif + +#ifdef DEBUG + write(iout,*) "------- SC restrs end -------" + write (iout,*) "------ After SC loop in e_modeller ------" + do i=loc_start,loc_end + write (iout,*) "i",i," gradc",(gradc(j,i,icg),j=1,3) + write (iout,*) "i",i," gradx",(gradx(j,i,icg),j=1,3) + enddo + if (waga_theta.eq.1.0d0) then + write (iout,*) "in e_modeller after SC restr end: dutheta" + do i=ithet_start,ithet_end + write (iout,*) i,dutheta(i) + enddo + endif + if (waga_d.eq.1.0d0) then + write (iout,*) "e_modeller after SC loop: duscdiff/x" + do i=1,nres + write (iout,*) i,(duscdiff(j,i),j=1,3) + write (iout,*) i,(duscdiffx(j,i),j=1,3) + enddo + endif +#endif + +! Total energy from homology restraints +#ifdef DEBUG + write (iout,*) "odleg",odleg," kat",kat +#endif +! +! Addition of energy of theta angle and SC local geom over constr_homologs ref strs +! +! ehomology_constr=odleg+kat +! +! For Lorentzian-type Urestr +! + + if (waga_dist.ge.0.0d0) then +! +! For Gaussian-type Urestr +! + ehomology_constr=(waga_dist*odleg+waga_angle*kat+ & + waga_theta*Eval+waga_d*Erot)*waga_homology(iset) +! write (iout,*) "ehomology_constr=",ehomology_constr +! print *, "ok7" + else +! +! For Lorentzian-type Urestr +! + ehomology_constr=(-waga_dist*odleg+waga_angle*kat+ & + waga_theta*Eval+waga_d*Erot)*waga_homology(iset) +! write (iout,*) "ehomology_constr=",ehomology_constr + print *, "ok8" + endif +#ifdef DEBUG + write (iout,*) "odleg",waga_dist,odleg," kat",waga_angle,kat, & + "Eval",waga_theta,eval, & + "Erot",waga_d,Erot + write (iout,*) "ehomology_constr",ehomology_constr +#endif + return +! +! FP 01/15 end +! + 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 subroutine e_modeller + +!---------------------------------------------------------------------------- subroutine ebend_kcc(etheta) logical lprn double precision thybt1(maxang_kcc),etheta @@ -11479,6 +12147,16 @@ enddo enddo +! write(iout,*), "const_homol",constr_homology + if (constr_homology.gt.0) then + do i=1,nct + do j=1,3 + gradc(j,i,icg)=gradc(j,i,icg)+duscdiff(j,i) +! write(iout,*) "duscdiff",duscdiff(j,i) + gradx(j,i,icg)=gradx(j,i,icg)+duscdiffx(j,i) + enddo + enddo + endif !#define DEBUG #ifdef DEBUG write (iout,*) "gloc before adding corr" @@ -12467,6 +13145,7 @@ subroutine check_ecartint ! Check the gradient of the energy in Cartesian coordinates. use io_base, only: intout + use MD_data, only: iset ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CONTROL' @@ -12499,6 +13178,7 @@ icg=1 nf=0 nfl=0 + if (iset.eq.0) iset=1 call intout ! call intcartderiv ! call checkintcartgrad @@ -12700,6 +13380,7 @@ subroutine check_ecartint ! Check the gradient of the energy in Cartesian coordinates. use io_base, only: intout + use MD_data, only: iset ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CONTROL' @@ -12732,6 +13413,7 @@ icg=1 nf=0 nfl=0 + if (iset.eq.0) iset=1 call intout ! call intcartderiv ! call checkintcartgrad @@ -12752,6 +13434,7 @@ enddo do j=1,3 grad_s(j,0)=gcart(j,0) + grad_s(j+3,0)=gxcart(j,0) enddo do i=1,nres do j=1,3 @@ -15973,7 +16656,7 @@ chip1=chip(itypi) integer :: i,n_corr,n_corr1,ierror,ierr real(kind=8) :: evdw2,evdw2_14,ehpb,etors,edihcnstr,etors_d,esccor,& evdw,ees,evdw1,eel_loc,eello_turn3,eello_turn4,& - ecorr,ecorr5,ecorr6,eturn6,time00 + ecorr,ecorr5,ecorr6,eturn6,time00, ehomology_constr ! write(iout,'(a,i2)')'Calling etotal_long ipot=',ipot !elwrite(iout,*)"in etotal long" @@ -15985,7 +16668,7 @@ chip1=chip(itypi) #endif endif !elwrite(iout,*)"in etotal long" - + ehomology_constr=0.0d0 #ifdef MPI ! write(iout,*) "ETOTAL_LONG Processor",fg_rank, ! & " absolute rank",myrank," nfgtasks",nfgtasks @@ -16183,6 +16866,7 @@ chip1=chip(itypi) energia(9)=eello_turn4 energia(10)=eturn6 energia(20)=Uconst+Uconst_back + energia(51)=ehomology_constr call sum_energy(energia,.true.) ! write (iout,*) "Exit ETOTAL_LONG" call flush(iout) @@ -16220,7 +16904,8 @@ chip1=chip(itypi) !el local variables integer :: i,nres6 real(kind=8) :: evdw,evdw1,evdw2,evdw2_14,esccor,etors_d,etors - real(kind=8) :: ehpb,escloc,estr,ebe,edihcnstr,ethetacnstr + real(kind=8) :: ehpb,escloc,estr,ebe,edihcnstr,ethetacnstr, & + ehomology_constr nres6=6*nres ! write(iout,'(a,i2)')'Calling etotal_short ipot=',ipot @@ -16436,6 +17121,16 @@ chip1=chip(itypi) call etor_d(etors_d) endif ! +! Homology restraints +! + if (constr_homology.ge.1) then + call e_modeller(ehomology_constr) +! print *,"tu" + else + ehomology_constr=0.0d0 + endif + +! ! 21/5/07 Calculate local sicdechain correlation energy ! if (wsccor.gt.0.0d0) then @@ -16470,6 +17165,7 @@ chip1=chip(itypi) energia(17)=estr energia(19)=edihcnstr energia(21)=esccor + energia(51)=ehomology_constr ! write (iout,*) "ETOTAL_SHORT before SUM_ENERGY" call flush(iout) call sum_energy(energia,.true.) @@ -16964,6 +17660,8 @@ chip1=chip(itypi) gvdwc_peppho(j,i)=0.0d0 gradnuclcatx(j,i)=0.0d0 gradnuclcat(j,i)=0.0d0 + duscdiff(j,i)=0.0d0 + duscdiffx(j,i)=0.0d0 enddo enddo do i=0,nres @@ -20201,7 +20899,7 @@ chip1=chip(itypi) else maxconts=10*nres ! (maxconts=maxres/4) endif - maxcont=12*nres ! Max. number of SC contacts + maxcont=100*nres ! Max. number of SC contacts maxvar=6*nres ! Max. number of variables !el maxdim=(nres-1)*(nres-2)/2 ! Max. number of derivatives of virtual-bond maxdim=nres*(nres-2)/2 ! Max. number of derivatives of virtual-bond @@ -20265,7 +20963,7 @@ chip1=chip(itypi) allocate(gacontm_hb3(3,maxconts,nres)) allocate(gacont_hbr(3,maxconts,nres)) allocate(grij_hb_cont(3,maxconts,nres)) -!(3,maxconts,maxres) + !(3,maxconts,maxres) allocate(facont_hb(maxconts,nres)) allocate(ees0p(maxconts,nres)) @@ -20539,8 +21237,8 @@ chip1=chip(itypi) allocate(dutheta(nres)) allocate(dugamma(nres)) !(maxres) - allocate(duscdiff(3,nres)) - allocate(duscdiffx(3,nres)) + allocate(duscdiff(3,-1:nres)) + allocate(duscdiffx(3,-1:nres)) !(3,maxres) !el i io:read_fragments ! allocate((:,:,:),allocatable :: wfrag_back !(3,maxfrag_back,maxprocs/20)