X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc%2Fenergy_p_new.F;h=652749c2b7abc4df9fc298025b797a77e82cddcb;hb=34d3ad3987785642be58fb2f26557d3314215577;hp=9b69cf70be79b01fd421d11892fcbe6d8a79bb83;hpb=f690e8b70bab14132839afebf080d4a28363b226;p=unres.git diff --git a/source/wham/src/energy_p_new.F b/source/wham/src/energy_p_new.F index 9b69cf7..652749c 100644 --- a/source/wham/src/energy_p_new.F +++ b/source/wham/src/energy_p_new.F @@ -107,7 +107,7 @@ c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t 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 @@ -116,7 +116,7 @@ c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t 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 @@ -154,6 +154,7 @@ c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t energia(19)=esccor energia(20)=edihcnstr energia(21)=evdw_t +c if (dyn_ss) call dyn_set_nss c detecting NaNQ #ifdef ISNAN #ifdef AIX @@ -770,6 +771,7 @@ C include 'COMMON.ENEPS' include 'COMMON.IOUNITS' include 'COMMON.CALC' + include 'COMMON.SBRIDGE' logical lprn common /srutu/icall integer icant @@ -800,6 +802,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 +883,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 +907,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 @@ -2869,24 +2888,16 @@ 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.ZSCOPT' include 'COMMON.SBRIDGE' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.VAR' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' - include 'COMMON.NAMES' dimension ggg(3) ehpb=0.0D0 -#ifdef DEBUG - do i=1,nres - write (iout,'(a4,2x,i4,3f10.5,5x,3f10.5)') restyp(itype(i)),i, - & (c(j,i),j=1,3),(c(j,i+nres),j=1,3) - enddo cd write(iout,*)'edis: nhpb=',nhpb,' fbr=',fbr cd write(iout,*)'link_start=',link_start,' link_end=',link_end -#endif if (link_end.eq.0) return do i=link_start,link_end C If ihpb(i) and jhpb(i) > NRES, this is a SC-SC distance, otherwise a @@ -2901,15 +2912,16 @@ C iii and jjj point to the residues for which the distance is assigned. iii=ii jjj=jj endif -#ifdef DEBUG - write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj, - & dhpb(i),dhpb1(i),forcon(i) -#endif +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 @@ -2917,10 +2929,8 @@ c Restraints from contact prediction 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 -#ifdef DEBUG - write (iout,*) "beta nmr", - & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) -#endif +c write (iout,*) "beta nmr", +c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) else dd=dist(ii,jj) rdis=dd-dhpb(i) @@ -2928,9 +2938,7 @@ C Get the force constant corresponding to this distance. waga=forcon(i) C Calculate the contribution to energy. ehpb=ehpb+waga*rdis*rdis -#ifdef DEBUG - write (iout,*) "beta reg",dd,waga*rdis*rdis -#endif +c write (iout,*) "beta reg",dd,waga*rdis*rdis C C Evaluate gradient. C @@ -2954,19 +2962,15 @@ C target distance. 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 -#ifdef DEBUG - write (iout,*) "alph nmr", - & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) -#endif +c write (iout,*) "alph nmr", +c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) else rdis=dd-dhpb(i) C Get the force constant corresponding to this distance. waga=forcon(i) C Calculate the contribution to energy. ehpb=ehpb+waga*rdis*rdis -#ifdef DEBUG - write (iout,*) "alpha reg",dd,waga*rdis*rdis -#endif +c write (iout,*) "alpha reg",dd,waga*rdis*rdis C C Evaluate gradient. C @@ -3050,11 +3054,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 @@ -3099,6 +3104,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 @@ -3118,8 +3124,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) @@ -3148,8 +3155,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) @@ -3429,6 +3437,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 @@ -3438,7 +3448,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 @@ -3452,13 +3462,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 @@ -3473,7 +3483,7 @@ c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1) enddo else phii1=0.0d0 - ityp3=nthetyp+1 + ityp3=ithetyp(itype(i)) do k=1,nsingle cosph2(k)=0.0d0 sinph2(k)=0.0d0 @@ -3582,10 +3592,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 @@ -4021,7 +4034,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 @@ -4556,6 +4570,8 @@ 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) @@ -4589,6 +4605,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 @@ -6474,7 +6493,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