X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fenergy.F90;h=ddc9833cd43ff65c32bf3957f2752600f8acda9e;hb=bc23440fbe68672d430f71f22f46b11265f003db;hp=63a76fb11c330f5f472172d0987d47bf6cdc4482;hpb=f7883bb11fc7df1ed921187999982670d26e6fa5;p=unres4.git diff --git a/source/unres/energy.F90 b/source/unres/energy.F90 index 63a76fb..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 @@ -400,12 +400,14 @@ if (nfgtasks.gt.1) then call MPI_Bcast(itime_mat,1,MPI_INT,king,FG_COMM,IERROR) endif + if (nres_molec(1).gt.0) then if (mod(itime_mat,imatupdate).eq.0) call make_SCp_inter_list ! write (iout,*) "after make_SCp_inter_list" if (mod(itime_mat,imatupdate).eq.0) call make_SCSC_inter_list ! write (iout,*) "after make_SCSC_inter_list" if (mod(itime_mat,imatupdate).eq.0) call make_pp_inter_list + endif ! write (iout,*) "after make_pp_inter_list" ! print *,'Processor',myrank,' calling etotal ipot=',ipot @@ -743,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 ! @@ -794,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 @@ -964,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" @@ -1005,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 @@ -1092,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 @@ -1105,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& @@ -1119,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& @@ -1257,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 @@ -1314,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,& @@ -1330,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)'/ & @@ -1378,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,& @@ -1387,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)'/ & @@ -1441,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 @@ -1915,7 +1932,7 @@ dGCLdOM1=0.0d0 dPOLdOM1=0.0d0 ! write (iout,*) "RWA", g_listscsc_start,g_listscsc_end,i,j - + if (nres_molec(1).eq.0) return do icont=g_listscsc_start,g_listscsc_end i=newcontlisti(icont) j=newcontlistj(icont) @@ -1951,7 +1968,7 @@ 'evdw',i,j,evdwij,' ss' ! if (energy_dec) write (iout,*) & ! 'evdw',i,j,evdwij,' ss' - do k=j+1,iend(i,iint) + do k=j+1,nres !C search over all next residues if (dyn_ss_mask(k)) then !C check if they are cysteins @@ -3448,6 +3465,7 @@ eel_loc=0.0d0 eello_turn3=0.0d0 eello_turn4=0.0d0 + if (nres_molec(1).eq.0) return ! if (icheckgrad.eq.1) then @@ -5061,14 +5079,19 @@ #ifdef NEWCORR gloc(nphi+i,icg)=gloc(nphi+i,icg)& -(gs13+gsE13+gsEE1)*wturn4& - *fac_shield(i)*fac_shield(j) + *fac_shield(i)*fac_shield(j) & + *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + gloc(nphi+i+1,icg)= gloc(nphi+i+1,icg)& -(gs23+gs21+gsEE2)*wturn4& - *fac_shield(i)*fac_shield(j) + *fac_shield(i)*fac_shield(j)& + *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) gloc(nphi+i+2,icg)= gloc(nphi+i+2,icg)& -(gs32+gsE31+gsEE3)*wturn4& - *fac_shield(i)*fac_shield(j) + *fac_shield(i)*fac_shield(j)& + *((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + !c gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)- !c & gs2 @@ -5394,6 +5417,7 @@ !d print '(a)','Enter ESCP' !d write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e ! do i=iatscp_s,iatscp_e + if (nres_molec(1).eq.0) return do icont=g_listscp_start,g_listscp_end i=newcontlistscpi(icont) j=newcontlistscpj(icont) @@ -7196,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) @@ -7542,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 @@ -11470,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" @@ -11761,7 +12448,7 @@ enddo do k=1,3 gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k)) -!C print *,'gg',k,gg(k) +! print *,'gg',k,gg(k) enddo ! print *,i,j,gg_lipi(3),gg_lipj(3),sss_ele_cut ! write (iout,*) "gg",(gg(k),k=1,3) @@ -12458,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' @@ -12490,6 +13178,7 @@ icg=1 nf=0 nfl=0 + if (iset.eq.0) iset=1 call intout ! call intcartderiv ! call checkintcartgrad @@ -12517,6 +13206,12 @@ do j=1,3 grad_s(j,i)=gcart(j,i) grad_s(j+3,i)=gxcart(j,i) + write(iout,*) "before movement analytical gradient" + do i=1,nres + write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),& + (gxcart(j,i),j=1,3) + enddo + enddo enddo else @@ -12685,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' @@ -12717,6 +13413,7 @@ icg=1 nf=0 nfl=0 + if (iset.eq.0) iset=1 call intout ! call intcartderiv ! call checkintcartgrad @@ -12737,16 +13434,20 @@ 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 grad_s(j,i)=gcart(j,i) -! if (i.eq.21) print *,"PRZEKAZANIE",gcart(j,i) - -! if (i.le.2) print *,"tu?!",gcart(j,i),grad_s(j,i),gxcart(j,i) grad_s(j+3,i)=gxcart(j,i) enddo enddo + write(iout,*) "before movement analytical gradient" + do i=1,nres + write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),& + (gxcart(j,i),j=1,3) + enddo + else !- split gradient check call zerograd @@ -15955,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" @@ -15967,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 @@ -16165,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) @@ -16202,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 @@ -16418,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 @@ -16452,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.) @@ -16712,13 +17426,13 @@ chip1=chip(itypi) #endif !#define DEBUG !el write (iout,*) "After sum_gradient" -#ifdef DEBUG - write (iout,*) "After sum_gradient" - do i=1,nres-1 - write (iout,*) i," gradc ",(gradc(j,i,icg),j=1,3) - write (iout,*) i," gradx ",(gradx(j,i,icg),j=1,3) - enddo -#endif +!#ifdef DEBUG +! write (iout,*) "After sum_gradient" +! do i=1,nres-1 +! write (iout,*) i," gradc ",(gradc(j,i,icg),j=1,3) +! write (iout,*) i," gradx ",(gradx(j,i,icg),j=1,3) +! enddo +!#endif !#undef DEBUG ! If performing constraint dynamics, add the gradients of the constraint energy if(usampl.and.totT.gt.eq_time) then @@ -16757,8 +17471,8 @@ chip1=chip(itypi) ! if (i.le.2) print *,"gcart_one",gcart(j,i),gradc(j,i,icg) enddo #ifdef DEBUG - write (iout,'(i5,2(3f10.5,5x),f10.5)') i,(gcart(j,i),j=1,3),& - (gxcart(j,i),j=1,3),gloc(i,icg) + write (iout,'(i5,2(3f10.5,5x),4f10.5)') i,(gcart(j,i),j=1,3),& + (gxcart(j,i),j=1,3),gloc(i,icg),(gloc_sc(j,i,icg),j=1,3) #endif enddo #ifdef TIMING @@ -16946,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 @@ -17087,10 +17803,12 @@ chip1=chip(itypi) do j=1,3 dcostheta(j,1,i)=-(dc_norm(j,i-1)+cost*dc_norm(j,i-2))/& vbld(i-1) - if (itype(i-1,1).ne.ntyp1) dtheta(j,1,i)=-dcostheta(j,1,i)/sint + if (((itype(i-1,1).ne.ntyp1).and.(sint.ne.0.0d0))) & + dtheta(j,1,i)=-dcostheta(j,1,i)/sint dcostheta(j,2,i)=-(dc_norm(j,i-2)+cost*dc_norm(j,i-1))/& vbld(i) - if (itype(i-1,1).ne.ntyp1) dtheta(j,2,i)=-dcostheta(j,2,i)/sint + if ((itype(i-1,1).ne.ntyp1).and.(sint.ne.0.0d0))& + dtheta(j,2,i)=-dcostheta(j,2,i)/sint enddo enddo #if defined(MPI) && defined(PARINTDER) @@ -17147,11 +17865,23 @@ chip1=chip(itypi) cost1=dcos(theta(i-1)) cosg=dcos(phi(i)) scalp=scalar(dc_norm(1,i-3),dc_norm(1,i-1)) + if ((sint*sint1).eq.0.0d0) then + fac0=0.0d0 + else fac0=1.0d0/(sint1*sint) + endif fac1=cost*fac0 fac2=cost1*fac0 + if (sint1.ne.0.0d0) then fac3=cosg*cost1/(sint1*sint1) + else + fac3=0.0d0 + endif + if (sint.ne.0.0d0) then fac4=cosg*cost/(sint*sint) + else + fac4=0.0d0 + endif ! Obtaining the gamma derivatives from sine derivative if (phi(i).gt.-pi4.and.phi(i).le.pi4.or. & phi(i).gt.pi34.and.phi(i).le.pi.or. & @@ -17160,8 +17890,16 @@ chip1=chip(itypi) call vecpr(dc_norm(1,i-3),dc_norm(1,i-1),vp2) call vecpr(dc_norm(1,i-3),dc_norm(1,i-2),vp3) do j=1,3 + if (sint.ne.0.0d0) then ctgt=cost/sint + else + ctgt=0.0d0 + endif + if (sint1.ne.0.0d0) then ctgt1=cost1/sint1 + else + ctgt1=0.0d0 + endif cosg_inv=1.0d0/cosg if (itype(i-1,1).ne.ntyp1 .and. itype(i-2,1).ne.ntyp1) then dsinphi(j,1,i)=-sing*ctgt1*dtheta(j,1,i-1) & @@ -17176,6 +17914,10 @@ chip1=chip(itypi) ! & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1) dphi(j,3,i)=cosg_inv*dsinphi(j,3,i) endif +! write(iout,*) "just after,close to pi",dphi(j,3,i),& +! sing*(ctgt1*dtheta(j,2,i-1)),ctgt*dtheta(j,1,i), & +! (fac0*vp2(j)+sing*dc_norm(j,i-2)),vbld_inv(i-1) + ! Bug fixed 3/24/05 (AL) enddo ! Obtaining the gamma derivatives from cosine derivative @@ -17230,12 +17972,25 @@ chip1=chip(itypi) ! write(iout,*) dc_norm2(j,i-2+nres),"dcnorm" enddo scalp=scalar(dc_norm2(1,i-2+nres),dc_norm(1,i-1)) + ! write(iout,*) "faki",fac0,fac1,fac2,fac3,fac + if ((sint*sint1).eq.0.0d0) then + fac0=0.0d0 + else fac0=1.0d0/(sint1*sint) + endif fac1=cost*fac0 fac2=cost1*fac0 + if (sint1.ne.0.0d0) then fac3=cosg*cost1/(sint1*sint1) + else + fac3=0.0d0 + endif + if (sint.ne.0.0d0) then fac4=cosg*cost/(sint*sint) - ! write(iout,*) "faki",fac0,fac1,fac2,fac3,fac4 + else + fac4=0.0d0 + endif + ! Obtaining the gamma derivatives from sine derivative if (tauangle(1,i).gt.-pi4.and.tauangle(1,i).le.pi4.or. & tauangle(1,i).gt.pi34.and.tauangle(1,i).le.pi.or. & @@ -17303,11 +18058,23 @@ chip1=chip(itypi) ! dc_norm2(j,i-1+nres)=-dc_norm(j,i-1+nres) ! enddo scalp=scalar(dc_norm(1,i-3),dc_norm(1,i-1+nres)) + if ((sint*sint1).eq.0.0d0) then + fac0=0.0d0 + else fac0=1.0d0/(sint1*sint) + endif fac1=cost*fac0 fac2=cost1*fac0 + if (sint1.ne.0.0d0) then fac3=cosg*cost1/(sint1*sint1) + else + fac3=0.0d0 + endif + if (sint.ne.0.0d0) then fac4=cosg*cost/(sint*sint) + else + fac4=0.0d0 + endif ! Obtaining the gamma derivatives from sine derivative if (tauangle(2,i).gt.-pi4.and.tauangle(2,i).le.pi4.or. & tauangle(2,i).gt.pi34.and.tauangle(2,i).le.pi.or. & @@ -17378,11 +18145,23 @@ chip1=chip(itypi) ! dc_norm2(j,i-1+nres)=-dc_norm(j,i-1+nres) enddo scalp=scalar(dc_norm2(1,i-2+nres),dc_norm(1,i-1+nres)) + if ((sint*sint1).eq.0.0d0) then + fac0=0.0d0 + else fac0=1.0d0/(sint1*sint) + endif fac1=cost*fac0 fac2=cost1*fac0 + if (sint1.ne.0.0d0) then fac3=cosg*cost1/(sint1*sint1) + else + fac3=0.0d0 + endif + if (sint.ne.0.0d0) then fac4=cosg*cost/(sint*sint) + else + fac4=0.0d0 + endif ! Obtaining the gamma derivatives from sine derivative if (tauangle(3,i).gt.-pi4.and.tauangle(3,i).le.pi4.or. & tauangle(3,i).gt.pi34.and.tauangle(3,i).le.pi.or. & @@ -19040,14 +19819,12 @@ chip1=chip(itypi) if (idssb(i).eq.newihpb(j) .and. & jdssb(i).eq.newjhpb(j)) found=.true. enddo -#ifndef CLUST -#ifndef WHAM +#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER) ! write(iout,*) "found",found,i,j if (.not.found.and.fg_rank.eq.0) & write(iout,'(a15,f12.2,f8.1,2i5)') & "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i) #endif -#endif enddo do i=1,newnss @@ -19057,21 +19834,22 @@ chip1=chip(itypi) if (newihpb(i).eq.idssb(j) .and. & newjhpb(i).eq.jdssb(j)) found=.true. enddo -#ifndef CLUST -#ifndef WHAM +#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER) ! write(iout,*) "found",found,i,j if (.not.found.and.fg_rank.eq.0) & write(iout,'(a15,f12.2,f8.1,2i5)') & "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i) #endif -#endif enddo - +!#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER) nss=newnss do i=1,nss idssb(i)=newihpb(i) jdssb(i)=newjhpb(i) enddo +!#else +! nss=0 +!#endif return end subroutine dyn_set_nss @@ -20121,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 @@ -20185,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)) @@ -20459,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) @@ -20551,10 +21329,10 @@ chip1=chip(itypi) !(3,3,2,maxres) ! allocateion of lists JPRDLA allocate(newcontlistppi(300*nres)) - allocate(newcontlistscpi(300*nres)) + allocate(newcontlistscpi(350*nres)) allocate(newcontlisti(300*nres)) allocate(newcontlistppj(300*nres)) - allocate(newcontlistscpj(300*nres)) + allocate(newcontlistscpj(350*nres)) allocate(newcontlistj(300*nres)) return @@ -22367,6 +23145,11 @@ chip1=chip(itypi) b3cav = alphasurcat(3,itypi,itypj) b4cav = alphasurcat(4,itypi,itypj) +! b1cav=0.0d0 +! b2cav=0.0d0 +! b3cav=0.0d0 +! b4cav=0.0d0 + ! used to determine whether we want to do quadrupole calculations eps_in = epsintabcat(itypi,itypj) if (eps_in.eq.0.0) eps_in=1.0 @@ -22402,7 +23185,7 @@ chip1=chip(itypi) enddo call to_box(chead(1,1),chead(2,1),chead(3,1)) call to_box(chead(1,2),chead(2,2),chead(3,2)) - +! write(iout,*) "TEST",chead(1,1),chead(2,1),chead(3,1),dc_norm(k, i+nres),d1 ! distance ! Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres)) ! Rsc(k) = Rsc_distance(k) * Rsc_distance(k) @@ -22455,6 +23238,16 @@ chip1=chip(itypi) rij_shift = Rtail - sig + sig0ij IF (rij_shift.le.0.0D0) THEN evdw = 1.0D20 + if (evdw.gt.1.0d6) then + write (*,'(2(1x,a3,i3),7f7.2)') & + restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,& + 1.0d0/rij,Rtail,Rhead,rij_shift, sig, sig0ij,sigsq + write(*,*) facsig,faceps1_inv,om1,chiom1,chi1 + write(*,*) "ANISO?!",chi1 +!evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,& +! Equad,evdwij+Fcav+eheadtail,evdw + endif + RETURN END IF sigder = -sig * sigsq @@ -22488,7 +23281,7 @@ chip1=chip(itypi) gg(1) = fac gg(2) = fac gg(3) = fac - +! print *,"GG(1),distance grad",gg(1) fac = chis1 * sqom1 + chis2 * sqom2 & - 2.0d0 * chis12 * om1 * om2 * om12 pom = 1.0d0 - chis1 * chis2 * sqom12 @@ -22527,12 +23320,12 @@ chip1=chip(itypi) erdxi = scalar( ertail(1), dC_norm(1,i+nres) ) erdxj = scalar( ertail(1), dC_norm(1,j) ) facd1 = dtailcat(1,itypi,itypj) * vbld_inv(i+nres) - facd2 = dtailcat(2,itypi,itypj) * vbld_inv(j+nres) + facd2 = dtailcat(2,itypi,itypj) * vbld_inv(j) DO k = 1, 3 pom = ertail(k)-facd1*(ertail(k)-erdxi*dC_norm(k,i+nres)) gradpepcatx(k,i) = gradpepcatx(k,i) & - (( dFdR + gg(k) ) * pom) - pom = ertail(k)-facd2*(ertail(k)-erdxj*dC_norm(k,j+nres)) + pom = ertail(k)-facd2*(ertail(k)-erdxj*dC_norm(k,j)) ! gvdwx(k,j) = gvdwx(k,j) & ! + (( dFdR + gg(k) ) * pom) gradpepcat(k,i) = gradpepcat(k,i) & @@ -22542,7 +23335,10 @@ chip1=chip(itypi) gg(k) = 0.0d0 ENDDO !c! Compute head-head and head-tail energies for each state +!! if (.false.) then ! turn off electrostatic + if (itype(j,5).gt.0) then ! the normal cation case isel = iabs(Qi) + 1 ! ion is always charged so iabs(Qj) +! print *,i,itype(i,1),isel IF (isel.eq.0) THEN !c! No charges - do nothing eheadtail = 0.0d0 @@ -22553,10 +23349,6 @@ chip1=chip(itypi) Qi=Qi*2 Qij=Qij*2 endif - if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then - Qj=Qj*2 - Qij=Qij*2 - endif CALL enq_cat(epol) eheadtail = epol @@ -22568,10 +23360,6 @@ chip1=chip(itypi) Qi=Qi*2 Qij=Qij*2 endif - if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then - Qj=Qj*2 - Qij=Qij*2 - endif ! write(iout,*) "KURWA0",d1 CALL edq_cat(ecl, elj, epol) @@ -22585,10 +23373,6 @@ chip1=chip(itypi) Qi=Qi*2 Qij=Qij*2 endif - if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then - Qj=Qj*2 - Qij=Qij*2 - endif CALL eqq_cat(Ecl,Egb,Epol,Fisocav,Elj) eheadtail = ECL + Egb + Epol + Fisocav + Elj @@ -22609,16 +23393,28 @@ chip1=chip(itypi) ! ! CALL energy_quad(istate,eheadtail,Ecl,Egb,Epol,Fisocav,Elj,Equad) END IF ! this endif ends the "catch the gly-gly" at the beggining of Fcav + else + write(iout,*) "not yet implemented",j,itype(j,5) + endif +!! endif ! turn off electrostatic evdw = evdw + Fcav + eheadtail +! if (evdw.gt.1.0d6) then +! write (*,'(2(1x,a3,i3),3f6.2,10f16.7)') & +! restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,& +! 1.0d0/rij,Rtail,Rhead,evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,& +! Equad,evdwij+Fcav+eheadtail,evdw +! endif IF (energy_dec) write (iout,'(2(1x,a3,i3),3f6.2,10f16.7)') & restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,& 1.0d0/rij,Rtail,Rhead,evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,& Equad,evdwij+Fcav+eheadtail,evdw ! evdw = evdw + Fcav + eheadtail - +! print *,"before sc_grad_cat", i,j, gradpepcat(1,j) ! iF (nstate(itypi,itypj).eq.1) THEN CALL sc_grad_cat +! print *,"after sc_grad_cat", i,j, gradpepcat(1,j) + ! END IF !c!------------------------------------------------------------------- !c! NAPISY KONCOWE @@ -22629,6 +23425,7 @@ chip1=chip(itypi) ! print *,"EVDW KURW",evdw,nres !!! return 17 continue +! go to 23 do i=ibond_start,ibond_end ! print *,"I am in EVDW",i @@ -22790,6 +23587,13 @@ chip1=chip(itypi) rij_shift = Rtail - sig + sig0ij IF (rij_shift.le.0.0D0) THEN evdw = 1.0D20 +! if (evdw.gt.1.0d6) then +! write (*,'(2(1x,a3,i3),6f6.2)') & +! restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,& +! 1.0d0/rij,Rtail,Rhead,rij_shift, sig, sig0ij +!evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,& +! Equad,evdwij+Fcav+eheadtail,evdw +! endif RETURN END IF sigder = -sig * sigsq @@ -22881,24 +23685,26 @@ chip1=chip(itypi) + (( dFdR + gg(k) ) * ertail(k)) gg(k) = 0.0d0 ENDDO + if (itype(j,5).gt.0) then !c! Compute head-head and head-tail energies for each state isel = 3 !c! Dipole-charge interactions - if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then - Qi=Qi*2 - Qij=Qij*2 - endif - if ((itype(j,1).eq.27).or.(itype(j,1).eq.26).or.(itype(j,1).eq.25)) then - Qj=Qj*2 - Qij=Qij*2 - endif CALL edq_cat_pep(ecl, elj, epol) eheadtail = ECL + elj + epol ! print *,"i,",i,eheadtail ! eheadtail = 0.0d0 - + else +!HERE WATER and other types of molecules solvents will be added + write(iout,*) "not yet implemented" +! CALL edd_cat_pep + endif evdw = evdw + Fcav + eheadtail - +! if (evdw.gt.1.0d6) then +! write (*,'(2(1x,a3,i3),3f6.2,10f16.7)') & +! restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,& +! 1.0d0/rij,Rtail,Rhead,evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,& +! Equad,evdwij+Fcav+eheadtail,evdw +! endif IF (energy_dec) write (iout,'(2(1x,a3,i3),3f6.2,10f16.7)') & restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,& 1.0d0/rij,Rtail,Rhead,evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,& @@ -22915,7 +23721,8 @@ chip1=chip(itypi) !c write (iout,*) "Number of loop steps in EGB:",ind !c energy_dec=.false. ! print *,"EVDW KURW",evdw,nres - + 23 continue +! print *,"before leave sc_grad_cat", i,j, gradpepcat(1,nres-1) return end subroutine ecats_prot_amber @@ -23470,9 +24277,13 @@ chip1=chip(itypi) dEvan2Cm,cm1,cm,vcat,vsug,v1,v2,dx,vcm,dEdipCm,dEdipCalp, & dEvan1Calp,dEvan2Cat,dEvan2Calp,dEtotalCat,dEdipCat,dEvan1Cat,dcosdcat, & dcosdcalp,dcosdcm,dEgbdCat,dEgbdCalp,dEgbdCm,dEcavdCat,dEcavdCalp, & - dEcavdCm + dEcavdCm,boxik real(kind=8),dimension(14) :: vcatnuclprm ecation_nucl=0.0d0 + boxik(1)=boxxsize + boxik(2)=boxysize + boxik(3)=boxzsize + if (nres_molec(5).eq.0) return itmp=0 do i=1,4 @@ -23515,12 +24326,16 @@ chip1=chip(itypi) vsug(k)=c(k,i) vcat(k)=c(k,j) enddo + call to_box(vcm(1),vcm(2),vcm(3)) + call to_box(vsug(1),vsug(2),vsug(3)) + call to_box(vcat(1),vcat(2),vcat(3)) do k=1,3 - dx(k) = vcat(k)-vcm(k) - enddo - do k=1,3 +! dx(k) = vcat(k)-vcm(k) +! enddo + dx(k)=boxshift(vcat(k)-vcm(k),boxik(k)) +! do k=1,3 v1(k)=dc(k,i+nres) - v2(k)=(vcat(k)-vsug(k)) + v2(k)=boxshift(vcat(k)-vsug(k),boxik(k)) enddo v1m = sqrt(v1(1)**2+v1(2)**2+v1(3)**2) v1dpdx = v1(1)*dx(1)+v1(2)*dx(2)+v1(3)*dx(3) @@ -25916,7 +26731,7 @@ chip1=chip(itypi) ! integer :: k !c! Epol and Gpol analytical parameters alphapol1 = alphapolcat(itypi,itypj) - alphapol2 = alphapolcat(itypj,itypi) + alphapol2 = alphapolcat2(itypj,itypi) !c! Fisocav and Gisocav analytical parameters al1 = alphisocat(1,itypi,itypj) al2 = alphisocat(2,itypi,itypj) @@ -26569,7 +27384,7 @@ chip1=chip(itypi) use calc_data use comm_momo double precision facd3, adler,epol - alphapol2 = alphapolcat(itypj,itypi) + alphapol2 = alphapolcat(itypi,itypj) !c! R2 - distance between head of jth side chain and tail of ith sidechain R2 = 0.0d0 DO k = 1, 3 @@ -26867,7 +27682,7 @@ chip1=chip(itypi) use calc_data double precision facd3, adler,ecl,elj,epol - alphapol2 = alphapolcat(itypj,itypi) + alphapol2 = alphapolcat(itypi,itypj) w1 = wqdipcat(1,itypi,itypj) w2 = wqdipcat(2,itypi,itypj) pis = sig0headcat(itypi,itypj) @@ -26986,7 +27801,7 @@ chip1=chip(itypi) use calc_data double precision facd3, adler,ecl,elj,epol - alphapol2 = alphapolcat(itypj,itypi) + alphapol2 = alphapolcat(itypi,itypj) w1 = wqdipcat(1,itypi,itypj) w2 = wqdipcat(2,itypi,itypj) pis = sig0headcat(itypi,itypj) @@ -27332,9 +28147,9 @@ chip1=chip(itypi) alf1 = 0.0d0 alf2 = 0.0d0 alf12 = 0.0d0 - dxj = dc_norm( 1, nres+j ) - dyj = dc_norm( 2, nres+j ) - dzj = dc_norm( 3, nres+j ) + dxj = 0.0d0 !dc_norm( 1, nres+j ) + dyj = 0.0d0 !dc_norm( 2, nres+j ) + dzj = 0.0d0 !dc_norm( 3, nres+j ) !c! distance from center of chain(?) to polar/charged head d1 = dheadcat(1, 1, itypi, itypj) d2 = dheadcat(2, 1, itypi, itypj) @@ -27584,8 +28399,10 @@ chip1=chip(itypi) zi=c(3,nres+i) call to_box(xi,yi,zi) do iint=1,nint_gr(i) +! print *,"is it wrong", iint,i do j=istart(i,iint),iend(i,iint) itypj=iabs(itype(j,1)) + if (energy_dec) write(iout,*) "LISTA ZAKRES",istart(i,iint),iend(i,iint),iatsc_s,iatsc_e if (itypj.eq.ntyp1) cycle xj=c(1,nres+j) yj=c(2,nres+j) @@ -27672,7 +28489,7 @@ chip1=chip(itypi) include 'mpif.h' real*8 :: xi,yi,zi,xj,yj,zj,xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp real*8 :: dist_init, dist_temp,r_buff_list - integer:: contlistscpi(250*nres),contlistscpj(250*nres) + integer:: contlistscpi(350*nres),contlistscpj(350*nres) ! integer :: newcontlistscpi(200*nres),newcontlistscpj(200*nres) integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint,ilist_scp,g_ilist_scp integer displ(0:nprocs),i_ilist_scp(0:nprocs),ierr