Merge branch 'prerelease-3.2.1' into czarek
[unres.git] / source / wham / src / energy_p_new.F
index 5ddb6f9..13fe796 100644 (file)
@@ -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
@@ -2897,9 +2916,12 @@ 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
       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
@@ -3081,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
@@ -3100,7 +3124,8 @@ c
           nbi=nbondterm(iti)
           if (nbi.eq.1) then
             diff=vbld(i+nres)-vbldsc0(1,iti)
-            write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),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
@@ -3130,7 +3155,8 @@ c
               usum=usum+uprod1
               usumsqder=usumsqder+ud(j)*uprod2
             enddo
-            write (iout,*) i,iti,vbld(i+nres),(vbldsc0(j,iti),
+            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
@@ -3564,10 +3590,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
@@ -4003,7 +4032,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