X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc%2Fenergy_p_new.F;h=1ac587d5dd811ad1998d6cbdcc07d3b0f678e67e;hb=c358637b84ca9f64f98d5205fe35dfbc916cbe61;hp=27f9010e803ad201d0282589d15b2f2486ccd980;hpb=715dd743e8e9985863b30940b432a0e466620ea9;p=unres.git diff --git a/source/wham/src/energy_p_new.F b/source/wham/src/energy_p_new.F index 27f9010..1ac587d 100644 --- a/source/wham/src/energy_p_new.F +++ b/source/wham/src/energy_p_new.F @@ -2,6 +2,7 @@ implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'DIMENSIONS.ZSCOPT' + include 'DIMENSIONS.FREE' #ifndef ISNAN external proc_proc @@ -22,6 +23,7 @@ cMS$ATTRIBUTES C :: proc_proc include 'COMMON.INTERACT' include 'COMMON.SBRIDGE' include 'COMMON.CHAIN' + include 'COMMON.CONTROL' double precision fact(6) cd write(iout, '(a,i2)')'Calling etotal ipot=',ipot cd print *,'nnt=',nnt,' nct=',nct @@ -102,25 +104,50 @@ c print *,ecorr,ecorr5,ecorr6,eturn6 if (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) then call multibody_hb(ecorr,ecorr5,ecorr6,n_corr,n_corr1) endif + + +c write(iout,*) "TEST_ENE1 constr_homology=",constr_homology + if (constr_homology.ge.1) then + call e_modeller(ehomology_constr) + else + ehomology_constr=0.0d0 + endif + +c write(iout,*) "TEST_ENE1 ehomology_constr=",ehomology_constr + +C BARTEK for dfa test! + if (wdfa_dist.gt.0) call edfad(edfadis) +c write(iout,*)'edfad is finished!', wdfa_dist,edfadis + if (wdfa_tor.gt.0) call edfat(edfator) +c write(iout,*)'edfat is finished!', wdfa_tor,edfator + if (wdfa_nei.gt.0) call edfan(edfanei) +c write(iout,*)'edfan is finished!', wdfa_nei,edfanei + if (wdfa_beta.gt.0) call edfab(edfabet) +c write(iout,*)'edfab is finished!', wdfa_beta,edfabet + c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t #ifdef SPLITELE etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2+welec*fact(1)*ees & +wvdwpp*evdw1 & +wang*ebe+wtor*fact(1)*etors+wscloc*escloc - & +wstrain*ehpb+nss*ebr+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5 + & +wstrain*ehpb+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5 & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4 & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6 & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d - & +wbond*estr+wsccor*fact(1)*esccor + & +wbond*estr+wsccor*fact(1)*esccor!+ehomology_constr + & +wdfa_dist*edfadis+wdfa_tor*edfator+wdfa_nei*edfanei + & +wdfa_beta*edfabet #else etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2 & +welec*fact(1)*(ees+evdw1) & +wang*ebe+wtor*fact(1)*etors+wscloc*escloc - & +wstrain*ehpb+nss*ebr+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5 + & +wstrain*ehpb+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5 & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4 & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6 & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d - & +wbond*estr+wsccor*fact(1)*esccor + & +wbond*estr+wsccor*fact(1)*esccor!+ehomology_constr + & +wdfa_dist*edfadis+wdfa_tor*edfator+wdfa_nei*edfanei + & +wdfa_beta*edfabet #endif energia(0)=etot energia(1)=evdw @@ -154,6 +181,12 @@ c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t energia(19)=esccor energia(20)=edihcnstr energia(21)=evdw_t + energia(22)=ehomology_constr + energia(23)=edfadis + energia(24)=edfator + energia(25)=edfanei + energia(26)=edfabet +c if (dyn_ss) call dyn_set_nss c detecting NaNQ #ifdef ISNAN #ifdef AIX @@ -191,7 +224,11 @@ C & wcorr5*fact(4)*gradcorr5(j,i)+ & wcorr6*fact(5)*gradcorr6(j,i)+ & wturn6*fact(5)*gcorr6_turn(j,i)+ - & wsccor*fact(2)*gsccorc(j,i) + & wsccor*fact(2)*gsccorc(j,i)+ + & wdfa_dist*gdfad(j,i)+ + & wdfa_tor*gdfat(j,i)+ + & wdfa_nei*gdfan(j,i)+ + & wdfa_beta*gdfab(j,i) gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+ & wbond*gradbx(j,i)+ & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+ @@ -210,7 +247,11 @@ C & wcorr5*fact(4)*gradcorr5(j,i)+ & wcorr6*fact(5)*gradcorr6(j,i)+ & wturn6*fact(5)*gcorr6_turn(j,i)+ - & wsccor*fact(2)*gsccorc(j,i) + & wsccor*fact(2)*gsccorc(j,i)+ + & wdfa_dist*gdfad(j,i)+ + & wdfa_tor*gdfat(j,i)+ + & wdfa_nei*gdfan(j,i)+ + & wdfa_beta*gdfab(j,i) gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+ & wbond*gradbx(j,i)+ & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+ @@ -268,6 +309,11 @@ C------------------------------------------------------------------------ esccor=energia(19) edihcnstr=energia(20) estr=energia(18) + ehomology_constr=energia(22) + edfadis=energia(23) + edfator=energia(24) + edfanei=energia(25) + edfabet=energia(26) #ifdef SPLITELE write (iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),evdw1, & wvdwpp, @@ -276,7 +322,9 @@ C------------------------------------------------------------------------ & ecorr,wcorr*fact(3),ecorr5,wcorr5*fact(4),ecorr6,wcorr6*fact(5), & eel_loc,wel_loc*fact(2),eello_turn3,wturn3*fact(2), & eello_turn4,wturn4*fact(3),eello_turn6,wturn6*fact(5), - & esccor,wsccor*fact(1),edihcnstr,ebr*nss,etot + & esccor,wsccor*fact(1),edihcnstr,ehomology_constr,ebr*nss, + & 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)'/ @@ -298,7 +346,12 @@ 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)'/ + & '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)') #else write (iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),estr,wbond, @@ -307,7 +360,9 @@ C------------------------------------------------------------------------ & ecorr6,wcorr6*fact(5),eel_loc,wel_loc*fact(2), & eello_turn3,wturn3*fact(2),eello_turn4,wturn4*fact(3), & eello_turn6,wturn6*fact(5),esccor*fact(1),wsccor, - & edihcnstr,ebr*nss,etot + & edihcnstr,ehomology_constr,ebr*nss, + & 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)'/ @@ -328,7 +383,12 @@ 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)'/ + & '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 @@ -770,6 +830,7 @@ C include 'COMMON.ENEPS' include 'COMMON.IOUNITS' include 'COMMON.CALC' + include 'COMMON.SBRIDGE' logical lprn common /srutu/icall integer icant @@ -800,6 +861,21 @@ C Calculate SC interaction energy. C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) +C in case of diagnostics write (iout,*) "TU SZUKAJ",i,j,dyn_ss_mask(i),dyn_ss_mask(j) +C /06/28/2013 Adasko: In case of dyn_ss - dynamic disulfide bond +C formation no electrostatic interactions should be calculated. If it +C would be allowed NaN would appear + IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN +C /06/28/2013 Adasko: dyn_ss_mask is logical statement wheather this Cys +C residue can or cannot form disulfide bond. There is still bug allowing +C Cys...Cys...Cys bond formation + call dyn_ssbond_ene(i,j,evdwij) +C /06/28/2013 Adasko: dyn_ssbond_ene is dynamic SS bond foration energy +C function in ssMD.F + evdw=evdw+evdwij +c if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') +c & 'evdw',i,j,evdwij,' ss' + ELSE ind=ind+1 itypj=itype(j) dscj_inv=vbld_inv(j+nres) @@ -866,6 +942,7 @@ c--------------------------------------------------------------- c write (iout,*) "i",i," j",j," itypi",itypi," itypj",itypj, c & " ij",ij," eneps",aux*e1/dabs(eps(itypi,itypj)), c & aux*e2/eps(itypi,itypj) +c write (iout,'(a6,2i5,0pf7.3)') 'evdw',i,j,evdwij if (lprn) then sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) epsi=bb(itypi,itypj)**2/aa(itypi,itypj) @@ -889,6 +966,7 @@ C Calculate the radial part of the gradient C Calculate angular part of the gradient. call sc_grad endif + ENDIF ! dyn_ss enddo ! j enddo ! iint enddo ! i @@ -1743,6 +1821,7 @@ C implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'DIMENSIONS.ZSCOPT' + include 'DIMENSIONS.FREE' include 'COMMON.CONTROL' include 'COMMON.IOUNITS' include 'COMMON.GEO' @@ -2869,12 +2948,14 @@ C Evaluate bridge-strain energy and its gradient in virtual-bond and SC vectors. C implicit real*8 (a-h,o-z) include 'DIMENSIONS' + include 'DIMENSIONS.FREE' include 'COMMON.SBRIDGE' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.VAR' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' + include 'COMMON.CONTROL' dimension ggg(3) ehpb=0.0D0 cd write(iout,*)'edis: nhpb=',nhpb,' fbr=',fbr @@ -2897,13 +2978,22 @@ c write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj, c & dhpb(i),dhpb1(i),forcon(i) C 24/11/03 AL: SS bridges handled separately because of introducing a specific C distance and angle dependent SS bond potential. + if (.not.dyn_ss .and. i.le.nss) then +C 15/02/13 CC dynamic SSbond - additional check if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then call ssbond_ene(iii,jjj,eij) ehpb=ehpb+2*eij + endif cd write (iout,*) "eij",eij else if (ii.gt.nres .and. jj.gt.nres) then c Restraints from contact prediction dd=dist(ii,jj) + if (constr_dist.eq.11) then + ehpb=ehpb+fordepth(i)**4.0d0 + & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)) + fac=fordepth(i)**4.0d0 + & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd + else if (dhpb1(i).gt.0.0d0) then ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd @@ -2921,7 +3011,8 @@ C C Evaluate gradient. C fac=waga*rdis/dd - endif + endif !end dhpb1(i).gt.0 + endif !end const_dist=11 do j=1,3 ggg(j)=fac*(c(j,jj)-c(j,ii)) enddo @@ -2937,6 +3028,20 @@ C C Calculate the distance between the two points and its difference from the C target distance. dd=dist(ii,jj) +C write(iout,*) "after",dd + if (constr_dist.eq.11) then + ehpb=ehpb+fordepth(i)**4.0d0 + & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)) + fac=fordepth(i)**4.0d0 + & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd +C ehpb=ehpb+fordepth(i)**4*rlornmr1(dd,dhpb(i),dhpb1(i)) +C fac=fordepth(i)**4*rlornmr1prim(dd,dhpb(i),dhpb1(i))/dd +C print *,ehpb,"tu?" +C write(iout,*) ehpb,"btu?", +C & dd,dhpb(i),dhpb1(i),fordepth(i),forcon(i) +C write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj, +C & ehpb,fordepth(i),dd + else if (dhpb1(i).gt.0.0d0) then ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd @@ -2954,6 +3059,7 @@ C Evaluate gradient. C fac=waga*rdis/dd endif + endif cd print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd, cd & ' waga=',waga,' fac=',fac do j=1,3 @@ -2974,7 +3080,7 @@ C Cartesian gradient in the SC vectors (ghpbx). enddo endif enddo - ehpb=0.5D0*ehpb + if (constr_dist.ne.11) ehpb=0.5D0*ehpb return end C-------------------------------------------------------------------------- @@ -3032,11 +3138,12 @@ C deltat12=om2-om1+2.0d0 cosphi=om12-om1*om2 eij=akcm*deltad*deltad+akth*(deltat1*deltat1+deltat2*deltat2) - & +akct*deltad*deltat12 + & +akct*deltad*deltat12+ebr +c & +akct*deltad*deltat12 & +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi -c write(iout,*) i,j,"rij",rij,"d0cm",d0cm," akcm",akcm," akth",akth, -c & " akct",akct," deltad",deltad," deltat",deltat1,deltat2, -c & " deltat12",deltat12," eij",eij + write(iout,*) i,j,"rij",rij,"d0cm",d0cm," akcm",akcm," akth",akth, + & " akct",akct," deltad",deltad," deltat",deltat1,deltat2, + & " deltat12",deltat12," eij",eij,"ebr",ebr ed=2*akcm*deltad+akct*deltat12 pom1=akct*deltad pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi @@ -3063,6 +3170,643 @@ C return end C-------------------------------------------------------------------------- +c MODELLER restraint function + subroutine e_modeller(ehomology_constr) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'DIMENSIONS.ZSCOPT' + include 'DIMENSIONS.FREE' + 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) + +c +c FP - 30/10/2014 Temporary specifications for homology restraints +c + double precision utheta_i,gutheta_i,sum_gtheta,sum_sgtheta, + & sgtheta + double precision, dimension (maxres) :: guscdiff,usc_diff + double precision, dimension (max_template) :: + & gtheta,dscdiff,uscdiffk,guscdiff2,guscdiff3, + & theta_diff + + 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.CONTROL' + include 'COMMON.HOMRESTR' +c + include 'COMMON.SETUP' + include 'COMMON.NAMES' + + do i=1,max_template + 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 +c 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) +c write (iout,*) "dij(",i,j,") =",dij + nexl=0 + do k=1,constr_homology + if(.not.l_homo(k,ii)) then + nexl=nexl+1 + cycle + endif + distance(k)=odl(k,ii)-dij +c write (iout,*) "distance(",k,") =",distance(k) +c +c For Gaussian-type Urestr +c + distancek(k)=0.5d0*distance(k)**2*sigma_odl(k,ii) ! waga_dist rmvd from Gaussian argument +c write (iout,*) "sigma_odl(",k,ii,") =",sigma_odl(k,ii) +c write (iout,*) "distancek(",k,") =",distancek(k) +c distancek(k)=0.5d0*waga_dist*distance(k)**2*sigma_odl(k,ii) +c +c For Lorentzian-type Urestr +c + 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 + +c min_odl=minval(distancek) + 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 +c 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 +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)) + if(.not.l_homo(k,ii)) cycle + if (waga_dist.ge.0.0d0) then +c +c For Gaussian-type Urestr +c + godl(k)=dexp(-distancek(k)+min_odl) + odleg2=odleg2+godl(k) +c +c For Lorentzian-type Urestr +c + else + odleg2=odleg2+distancek(k) + endif + +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 +c write (iout,*) "godl",(godl(k),k=1,constr_homology) ! exponents +c 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 +c +c For Gaussian-type Urestr +c + odleg=odleg-dLOG(odleg2/constr_homology)+min_odl +c +c For Lorentzian-type Urestr +c + else + odleg=odleg+odleg2/constr_homology + endif +c +#ifdef GRAD +c write (iout,*) "odleg",odleg ! sum of -ln-s +c Gradient +c +c For Gaussian-type Urestr +c + if (waga_dist.ge.0.0d0) sum_godl=odleg2 + sum_sgodl=0.0d0 + do k=1,constr_homology +c godl=dexp(((-(distance(i,j,k)**2)/(2*(sigma_odl(i,j,k))**2)) +c & *waga_dist)+min_odl +c sgodl=-godl(k)*distance(k)*sigma_odl(k,ii)*waga_dist +c + if(.not.l_homo(k,ii)) cycle + if (waga_dist.ge.0.0d0) then +c For Gaussian-type Urestr +c + sgodl=-godl(k)*distance(k)*sigma_odl(k,ii) ! waga_dist rmvd +c +c For Lorentzian-type Urestr +c + else + sgodl=-2*sigma_odlir(k,ii)*(distance(k)/(distance(k)**2+ + & sigma_odlir(k,ii)**2)**2) + endif + 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 + if (waga_dist.ge.0.0d0) then +c +c For Gaussian-type Urestr +c + grad_odl3=waga_homology(iset)*waga_dist + & *sum_sgodl/(sum_godl*dij) +c +c For Lorentzian-type Urestr +c + else +c Original grad expr modified by analogy w Gaussian-type Urestr grad +c grad_odl3=-waga_homology(iset)*waga_dist*sum_sgodl + grad_odl3=-waga_homology(iset)*waga_dist* + & sum_sgodl/(constr_homology*dij) + endif +c +c 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) +c if (i.eq.25.and.j.eq.27) then +c write(iout,*) "jik",jik,"i",i,"j",j +c write(iout,*) "sum_sgodl",sum_sgodl,"sgodl",sgodl +c write(iout,*) "grad_odl3",grad_odl3 +c write(iout,*) "c(",jik,i,")",c(jik,i),"c(",jik,j,")",c(jik,j) +c write(iout,*) "ggodl",ggodl +c write(iout,*) "ghpbc(",jik,i,")", +c & ghpbc(jik,i),"ghpbc(",jik,j,")", +c & ghpbc(jik,j) +c endif + enddo +#endif +ccc write(iout,778)"TEST: odleg2=", odleg2, "DLOG(odleg2)=", +ccc & dLOG(odleg2),"-odleg=", -odleg + + enddo ! ii-loop for dist +#ifdef DEBUG + write(iout,*) "------- dist restrs end -------" +c if (waga_angle.eq.1.0d0 .or. waga_theta.eq.1.0d0 .or. +c & waga_d.eq.1.0d0) call sum_gradient +#endif +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 +#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 +#ifdef OLDRESTR + kat2=0.0d0 +#else + kat2=nexl +#endif +c betai=beta(i,i+1,i+2,i+3) + betai = phi(i) +c write (iout,*) "betai =",betai + do k=1,constr_homology + dih_diff(k)=pinorm(dih(k,i)-betai) +c write (iout,*) "dih_diff(",k,") =",dih_diff(k) +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) +#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) +#endif +c kat3=-0.5d0*waga_angle*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 +c write (iout,*) "gdih",(gdih(k),k=1,constr_homology) ! exps +c 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) +c write (iout,*) "kat",kat ! sum of -ln-s + +ccc write(iout,778)"TEST: kat2=", kat2, "DLOG(kat2)=", +ccc & dLOG(kat2), "-kat=", -kat + +#ifdef GRAD +c ---------------------------------------------------------------------- +c Gradient +c ---------------------------------------------------------------------- + + 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) +#endif +c sgdih=-gdih(k)*dih_diff(k)*sigma_dih(k,i)*waga_angle + sum_sgdih=sum_sgdih+sgdih + enddo +c grad_dih3=sum_sgdih/sum_gdih + grad_dih3=waga_homology(iset)*waga_angle*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 +c if (i.eq.25) then +c write(iout,*) "i",i,"icg",icg,"gloc(",i,icg,")",gloc(i,icg) +c endif +ccc write(iout,747) "GRAD_KAT_2", i, nphi, icg, grad_dih3, +ccc & gloc(nphi+i-3,icg) +#endif + enddo ! i-loop for dih +#ifdef DEBUG + write(iout,*) "------- dih restrs end -------" +#endif + +c Pseudo-energy and gradient for theta angle restraints from +c homology templates +c FP 01/15 - inserted from econstr_local_test.F, loop structure +c adapted + +c +c For constr_homology reference structures (FP) +c +c Uconst_back_tot=0.0d0 + Eval=0.0d0 + Erot=0.0d0 +c Econstr_back legacy +#ifdef GRAD + do i=1,nres +c do i=ithet_start,ithet_end + dutheta(i)=0.0d0 +c enddo +c do i=loc_start,loc_end + do j=1,3 + duscdiff(j,i)=0.0d0 + duscdiffx(j,i)=0.0d0 + enddo + enddo +#endif +c +c do iref=1,nref +c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end +c write (iout,*) "waga_theta",waga_theta + if (waga_theta.gt.0.0d0) then +#ifdef DEBUG + write (iout,*) "usampl",usampl + write(iout,*) "------- theta restrs start -------" +c do i=ithet_start,ithet_end +c write (iout,*) "gloc_init(",nphi+i,icg,")",gloc(nphi+i,icg) +c enddo +#endif +c write (iout,*) "maxres",maxres,"nres",nres + + do i=ithet_start,ithet_end +c +c do i=1,nfrag_back +c ii = ifrag_back(2,i,iset)-ifrag_back(1,i,iset) +c +c Deviation of theta angles wrt constr_homology ref structures +c + utheta_i=0.0d0 ! argument of Gaussian for single k +#ifdef OLDRESTR + gutheta_i=0.0d0 ! Sum of Gaussians over constr_homology ref structures +#else + gutheta_i=nexl +#endif +c do j=ifrag_back(1,i,iset)+2,ifrag_back(2,i,iset) ! original loop +c over residues in a fragment +c write (iout,*) "theta(",i,")=",theta(i) + do k=1,constr_homology +c +c dtheta_i=theta(j)-thetaref(j,iref) +c dtheta_i=thetaref(k,i)-theta(i) ! original form without indexing + theta_diff(k)=thetatpl(k,i)-theta(i) +c + utheta_i=-0.5d0*theta_diff(k)**2*sigma_theta(k,i) ! waga_theta rmvd from Gaussian argument +c 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+dexp(utheta_i) ! Sum of Gaussians (pk) +c Gradient for single Gaussian restraint in subr Econstr_back +c dutheta(j-2)=dutheta(j-2)+wfrag_back(1,i,iset)*dtheta_i/(ii-1) +c + enddo +c write (iout,*) "gtheta",(gtheta(k),k=1,constr_homology) ! exps +c write (iout,*) "i",i," gutheta_i",gutheta_i ! sum of exps + +c +#ifdef GRAD +c Gradient for multiple Gaussian restraint + sum_gtheta=gutheta_i + sum_sgtheta=0.0d0 + do k=1,constr_homology +c New generalized expr for multiple Gaussian from Econstr_back + sgtheta=-gtheta(k)*theta_diff(k)*sigma_theta(k,i) ! waga_theta rmvd +c +c sgtheta=-gtheta(k)*theta_diff(k)*sigma_theta(k,i)*waga_theta ! right functional form? + sum_sgtheta=sum_sgtheta+sgtheta ! cum variable + enddo +c +c Final value of gradient using same var as in Econstr_back + dutheta(i-2)=sum_sgtheta/sum_gtheta*waga_theta + & *waga_homology(iset) +c dutheta(i)=sum_sgtheta/sum_gtheta +c +c Uconst_back=Uconst_back+waga_theta*utheta(i) ! waga_theta added as weight +#endif + Eval=Eval-dLOG(gutheta_i/constr_homology) +c write (iout,*) "utheta(",i,")=",utheta(i) ! -ln of sum of exps +c write (iout,*) "Uconst_back",Uconst_back ! sum of -ln-s +c Uconst_back=Uconst_back+utheta(i) + enddo ! (i-loop for theta) +#ifdef DEBUG + write(iout,*) "------- theta restrs end -------" +#endif + endif +c +c Deviation of local SC geometry +c +c Separation of two i-loops (instructed by AL - 11/3/2014) +c +c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end +c 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 +#ifdef OLDRESTR + guscdiff(i)=0.0d0 ! Sum of Gaussians over constr_homology ref structures +#else + guscdiff(i)=nexl +#endif +c do j=ifrag_back(1,i,iset)+1,ifrag_back(2,i,iset)-1 ! Econstr_back legacy +c write(iout,*) "xxtab, yytab, zztab" +c write(iout,'(i5,3f8.2)') i,xxtab(i),yytab(i),zztab(i) + do k=1,constr_homology +c + dxx=-xxtpl(k,i)+xxtab(i) ! Diff b/w x component of ith SC vector in model and kth ref str? +c 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 +c write(iout,*) "dxx, dyy, dzz" +c write(iout,'(2i5,3f8.2)') k,i,dxx,dyy,dzz +c + usc_diff_i=-0.5d0*(dxx**2+dyy**2+dzz**2)*sigma_d(k,i) ! waga_d rmvd from Gaussian argument +c usc_diff(i)=-0.5d0*waga_d*(dxx**2+dyy**2+dzz**2)*sigma_d(k,i) ! waga_d? +c uscdiffk(k)=usc_diff(i) + guscdiff2(k)=dexp(usc_diff_i) ! without min_scdiff + guscdiff(i)=guscdiff(i)+dexp(usc_diff_i) !Sum of Gaussians (pk) +c write (iout,'(i5,6f10.5)') j,xxtab(j),yytab(j),zztab(j), +c & xxref(j),yyref(j),zzref(j) + enddo +c +c Gradient +c +c Generalized expression for multiple Gaussian acc to that for a single +c Gaussian in Econstr_back as instructed by AL (FP - 03/11/2014) +c +c Original implementation +c sum_guscdiff=guscdiff(i) +c +c sum_sguscdiff=0.0d0 +c do k=1,constr_homology +c sguscdiff=-guscdiff2(k)*dscdiff(k)*sigma_d(k,i)*waga_d !waga_d? +c sguscdiff=-guscdiff3(k)*dscdiff(k)*sigma_d(k,i)*waga_d ! w min_uscdiff +c sum_sguscdiff=sum_sguscdiff+sguscdiff +c enddo +c +c Implementation of new expressions for gradient (Jan. 2015) +c +c grad_uscdiff=sum_sguscdiff/(sum_guscdiff*dtab) !? +#ifdef GRAD + do k=1,constr_homology +c +c New calculation of dxx, dyy, and dzz corrected by AL (07/11), was missing and wrong +c before. Now the drivatives should be correct +c + dxx=-xxtpl(k,i)+xxtab(i) ! Diff b/w x component of ith SC vector in model and kth ref str? +c 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 +c +c New implementation +c + sum_guscdiff=guscdiff2(k)*!(dsqrt(dxx*dxx+dyy*dyy+dzz*dzz))* -> wrong! + & sigma_d(k,i) ! for the grad wrt r' +c sum_sguscdiff=sum_sguscdiff+sum_guscdiff +c +c +c 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) +c +#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) +c write(iout,*) "sum_sguscdiff",sum_sguscdiff +cc write(iout,*) "dXX_Ctab(",jik,i,")",dXX_Ctab(jik,i) +c write(iout,*) "dYY_Ctab(",jik,i,")",dYY_Ctab(jik,i) +c write(iout,*) "dZZ_Ctab(",jik,i,")",dZZ_Ctab(jik,i) +c write(iout,*) "dXX_C1tab(",jik,i,")",dXX_C1tab(jik,i) +c write(iout,*) "dYY_C1tab(",jik,i,")",dYY_C1tab(jik,i) +c write(iout,*) "dZZ_C1tab(",jik,i,")",dZZ_C1tab(jik,i) +c write(iout,*) "dXX_XYZtab(",jik,i,")",dXX_XYZtab(jik,i) +c write(iout,*) "dYY_XYZtab(",jik,i,")",dYY_XYZtab(jik,i) +c write(iout,*) "dZZ_XYZtab(",jik,i,")",dZZ_XYZtab(jik,i) +c write(iout,*) "duscdiff(",jik,i-1,")",duscdiff(jik,i-1) +c write(iout,*) "duscdiff(",jik,i,")",duscdiff(jik,i) +c write(iout,*) "duscdiffx(",jik,i,")",duscdiffx(jik,i) +c endif +#endif + enddo + enddo +#endif +c +c uscdiff(i)=-dLOG(guscdiff(i)/(ii-1)) ! Weighting by (ii-1) required? +c usc_diff(i)=-dLOG(guscdiff(i)/constr_homology) ! + min_uscdiff ? +c +c write (iout,*) i," uscdiff",uscdiff(i) +c +c Put together deviations from local geometry + +c Uconst_back=Uconst_back+wfrag_back(1,i,iset)*utheta(i)+ +c & wfrag_back(3,i,iset)*uscdiff(i) + Erot=Erot-dLOG(guscdiff(i)/constr_homology) +c write (iout,*) "usc_diff(",i,")=",usc_diff(i) ! -ln of sum of exps +c write (iout,*) "Uconst_back",Uconst_back ! cum sum of -ln-s +c Uconst_back=Uconst_back+usc_diff(i) +c +c Gradient of multiple Gaussian restraint (FP - 04/11/2014 - right?) +c +c New implment: multiplied by sum_sguscdiff +c + + enddo ! (i-loop for dscdiff) + +c 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 + +c Total energy from homology restraints +#ifdef DEBUG + write (iout,*) "odleg",odleg," kat",kat + write (iout,*) "odleg",odleg," kat",kat + write (iout,*) "Eval",Eval," Erot",Erot + write (iout,*) "waga_homology(",iset,")",waga_homology(iset) + write (iout,*) "waga_dist ",waga_dist,"waga_angle ",waga_angle + write (iout,*) "waga_theta ",waga_theta,"waga_d ",waga_d +#endif +c +c Addition of energy of theta angle and SC local geom over constr_homologs ref strs +c +c ehomology_constr=odleg+kat +c +c For Lorentzian-type Urestr +c + + if (waga_dist.ge.0.0d0) then +c +c For Gaussian-type Urestr +c +c ehomology_constr=(waga_dist*odleg+waga_angle*kat+ +c & waga_theta*Eval+waga_d*Erot)*waga_homology(iset) + ehomology_constr=waga_dist*odleg+waga_angle*kat+ + & waga_theta*Eval+waga_d*Erot +c write (iout,*) "ehomology_constr=",ehomology_constr + else +c +c For Lorentzian-type Urestr +c +c ehomology_constr=(-waga_dist*odleg+waga_angle*kat+ +c & waga_theta*Eval+waga_d*Erot)*waga_homology(iset) + ehomology_constr=-waga_dist*odleg+waga_angle*kat+ + & waga_theta*Eval+waga_d*Erot +c write (iout,*) "ehomology_constr=",ehomology_constr + 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 + + 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 ebond(estr) c c Evaluate the energy of stretching of the CA-CA and CA-SC virtual bonds @@ -3070,6 +3814,7 @@ c implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'DIMENSIONS.ZSCOPT' + include 'DIMENSIONS.FREE' include 'COMMON.LOCAL' include 'COMMON.GEO' include 'COMMON.INTERACT' @@ -3081,6 +3826,7 @@ c include 'COMMON.FFIELD' include 'COMMON.CONTROL' double precision u(3),ud(3) + logical :: lprn=.false. estr=0.0d0 do i=nnt+1,nct diff = vbld(i)-vbldp0 @@ -3100,8 +3846,9 @@ c nbi=nbondterm(iti) if (nbi.eq.1) then diff=vbld(i+nres)-vbldsc0(1,iti) -c write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),diff, -c & AKSC(1,iti),AKSC(1,iti)*diff*diff + if (lprn) + & write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),diff, + & AKSC(1,iti),AKSC(1,iti)*diff*diff estr=estr+0.5d0*AKSC(1,iti)*diff*diff do j=1,3 gradbx(j,i)=AKSC(1,iti)*diff*dc(j,i+nres)/vbld(i+nres) @@ -3130,8 +3877,9 @@ c & AKSC(1,iti),AKSC(1,iti)*diff*diff usum=usum+uprod1 usumsqder=usumsqder+ud(j)*uprod2 enddo -c write (iout,*) i,iti,vbld(i+nres),(vbldsc0(j,iti), -c & AKSC(j,iti),abond0(j,iti),u(j),j=1,nbi) + if (lprn) + & write (iout,*) i,iti,vbld(i+nres),(vbldsc0(j,iti), + & AKSC(j,iti),abond0(j,iti),u(j),j=1,nbi) estr=estr+uprod/usum do j=1,3 gradbx(j,i)=usumsqder/(usum*usum)*dc(j,i+nres)/vbld(i+nres) @@ -3393,6 +4141,7 @@ C implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'DIMENSIONS.ZSCOPT' + include 'DIMENSIONS.FREE' include 'COMMON.LOCAL' include 'COMMON.GEO' include 'COMMON.INTERACT' @@ -3411,6 +4160,8 @@ C etheta=0.0D0 c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1) do i=ithet_start,ithet_end + if ((itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or. + & (itype(i).eq.ntyp1)) cycle dethetai=0.0d0 dephii=0.0d0 dephii1=0.0d0 @@ -3420,7 +4171,7 @@ c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1) coskt(k)=dcos(k*theti2) sinkt(k)=dsin(k*theti2) enddo - if (i.gt.3) then + if (i.gt.3 .and. itype(max0(i-3,1)).ne.ntyp1) then #ifdef OSF phii=phi(i) if (phii.ne.phii) phii=150.0 @@ -3434,13 +4185,13 @@ c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1) enddo else phii=0.0d0 - ityp1=nthetyp+1 + ityp1=ithetyp(itype(i-2)) do k=1,nsingle cosph1(k)=0.0d0 sinph1(k)=0.0d0 enddo endif - if (i.lt.nres) then + if (i.lt.nres .and. itype(i+1).ne.ntyp1) then #ifdef OSF phii1=phi(i+1) if (phii1.ne.phii1) phii1=150.0 @@ -3455,7 +4206,8 @@ c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1) enddo else phii1=0.0d0 - ityp3=nthetyp+1 +c ityp3=nthetyp+1 + ityp3=ithetyp(itype(i)) do k=1,nsingle cosph2(k)=0.0d0 sinph2(k)=0.0d0 @@ -3564,10 +4316,13 @@ c call flush(iout) 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 @@ -3868,6 +4623,7 @@ C implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'DIMENSIONS.ZSCOPT' + include 'DIMENSIONS.FREE' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.VAR' @@ -4003,7 +4759,8 @@ c sumene = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1)) escloc = escloc + sumene c write (2,*) "escloc",escloc if (.not. calc_grad) goto 1 -#ifdef DEBUG + +#ifdef DEBUG2 C C This section to check the numerical derivatives of the energy of ith side C chain in xx, yy, zz, and theta. Use the -DDEBUG compiler option or insert @@ -4518,6 +5275,7 @@ c amino-acid residues. implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'DIMENSIONS.ZSCOPT' + include 'DIMENSIONS.FREE' include 'COMMON.VAR' include 'COMMON.GEO' include 'COMMON.LOCAL' @@ -4534,10 +5292,11 @@ c amino-acid residues. C Set lprn=.true. for debugging lprn=.false. c lprn=.true. -c write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor +c write (iout,*) "EBACK_SC_COR",itau_start,itau_end,nterm_sccor esccor=0.0D0 do i=itau_start,itau_end esccor_ii=0.0D0 + if ((itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)) cycle isccori=isccortyp(itype(i-2)) isccori1=isccortyp(itype(i-1)) phii=phi(i) @@ -4571,6 +5330,9 @@ c 3 = SC...Ca...Ca...SCi cosphi=dcos(j*tauangle(intertyp,i)) sinphi=dsin(j*tauangle(intertyp,i)) esccor=esccor+v1ij*cosphi+v2ij*sinphi +#ifdef DEBUG + esccor_ii=esccor_ii+v1ij*cosphi+v2ij*sinphi +#endif gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi) enddo gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci @@ -4583,6 +5345,9 @@ c &gloc_sc(intertyp,i-3,icg) & ,(v2sccor(j,intertyp,itori,itori1),j=1,6) gsccor_loc(i-3)=gsccor_loc(i-3)+gloci enddo !intertyp +#ifdef DEBUG + write (iout,*) "i",i,(tauangle(j,i),j=1,3),esccor_ii +#endif enddo c do i=1,nres c write (iout,*) "W@T@F", gloc_sc(1,i,icg),gloc(i,icg) @@ -4904,6 +5669,10 @@ c & ' jj=',jj,' kk=',kk C Contacts I-J and (I+1)-(J+1) or (I+1)-(J-1) occur simultaneously. C The system gains extra energy. ecorr=ecorr+ehbcorr(i,j,i+1,j1,jj,kk,0.72D0,0.32D0) +#ifdef DEBUG + write (iout,*) "ecorr",i,j,i+1,j1, + & ehbcorr(i,j,i+1,j1,jj,kk,0.72D0,0.32D0) +#endif n_corr=n_corr+1 else if (j1.eq.j) then C Contacts I-J and I-(J+1) occur simultaneously. @@ -5191,11 +5960,11 @@ cd ees0pkl=0.0D0 cd ees0pij=1.0D0 cd ees0mkl=0.0D0 cd ees0mij=1.0D0 -c write (iout,*)'Contacts have occurred for peptide groups',i,j, -c & ' and',k,l -c write (iout,*)'Contacts have occurred for peptide groups', -c & i,j,' fcont:',eij,' eij',' eesij',ees0pij,ees0mij,' and ',k,l -c & ,' fcont ',ekl,' eeskl',ees0pkl,ees0mkl,' ees=',ees +cd write (iout,*)'Contacts have occurred for peptide groups',i,j, +cd & ' and',k,l +cd write (iout,*)'Contacts have occurred for peptide groups', +cd & i,j,' fcont:',eij,' eij',' eesij',ees0pij,ees0mij,' and ',k,l +cd & ,' fcont ',ekl,' eeskl',ees0pkl,ees0mkl,' ees=',ees C Calculate the multi-body contribution to energy. ecorr=ecorr+ekont*ees if (calc_grad) then @@ -6456,7 +7225,7 @@ c---------------------------------------------------------------------------- include 'COMMON.GEO' logical swap double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2), - & auxvec1(2),auxvec2(1),auxmat1(2,2) + & auxvec1(2),auxvec2(2),auxmat1(2,2) logical lprn common /kutas/ lprn CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC