update
[unres.git] / source / cluster / wham / src / energy_p_new.F
index ee00811..1c4ef0c 100644 (file)
@@ -107,7 +107,7 @@ C
 #ifdef SPLITELE
       etot=wsc*evdw+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
@@ -115,7 +115,7 @@ C
 #else
       etot=wsc*evdw+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
@@ -152,6 +152,7 @@ C
       energia(18)=estr
       energia(19)=esccor
       energia(20)=edihcnstr
+cc      if (dyn_ss) call dyn_set_nss
 c detecting NaNQ
       i=0
 #ifdef WINPGI
@@ -218,7 +219,7 @@ cd        write (iout,*) i,g_corr5_loc(i)
      &   +wturn4*fact(3)*gel_loc_turn4(i)
      &   +wturn3*fact(2)*gel_loc_turn3(i)
      &   +wturn6*fact(5)*gel_loc_turn6(i)
-     &   +wel_loc*fact(2)*gel_loc_loc(i)+
+     &   +wel_loc*fact(2)*gel_loc_loc(i)
      &   +wsccor*fact(1)*gsccor_loc(i)
       enddo
       endif
@@ -723,6 +724,7 @@ c      include "DIMENSIONS.COMPAR"
       include 'COMMON.INTERACT'
       include 'COMMON.IOUNITS'
       include 'COMMON.CALC'
+      include 'COMMON.SBRIDGE'
       logical lprn
       common /srutu/icall
       integer icant
@@ -748,6 +750,12 @@ C Calculate SC interaction energy.
 C
         do iint=1,nint_gr(i)
           do j=istart(i,iint),iend(i,iint)
+            IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
+              call dyn_ssbond_ene(i,j,evdwij)
+              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)
@@ -830,6 +838,7 @@ C Calculate the radial part of the gradient
 C Calculate angular part of the gradient.
             call sc_grad
             endif
+            ENDIF    ! SSBOND
           enddo      ! j
         enddo        ! iint
       enddo          ! i
@@ -854,6 +863,7 @@ c      include "DIMENSIONS.COMPAR"
       include 'COMMON.INTERACT'
       include 'COMMON.IOUNITS'
       include 'COMMON.CALC'
+      include 'COMMON.SBRIDGE'
       common /srutu/ icall
       logical lprn
       integer icant
@@ -879,6 +889,13 @@ 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)
+            IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
+              call dyn_ssbond_ene(i,j,evdwij)
+              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)
@@ -961,6 +978,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
       include 'COMMON.VAR'
       include 'COMMON.INTERACT'
       include 'COMMON.IOUNITS'
-      include 'COMMON.NAMES'
       dimension ggg(3)
       ehpb=0.0D0
 cd      write(iout,*)'edis: nhpb=',nhpb,' fbr=',fbr
 cd      write(iout,*)'link_start=',link_start,' link_end=',link_end
-#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
-#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
@@ -2825,26 +2836,25 @@ 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
 cd          write (iout,*) "eij",eij
+        endif
         else if (ii.gt.nres .and. jj.gt.nres) then
 c Restraints from contact prediction
           dd=dist(ii,jj)
           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)
@@ -2852,9 +2862,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
@@ -2878,19 +2886,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
@@ -2974,7 +2978,7 @@ 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
      &  +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,
@@ -3351,6 +3355,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
@@ -3360,7 +3366,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
@@ -3374,13 +3380,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
@@ -3395,7 +3401,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
@@ -4505,6 +4511,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
@@ -4886,6 +4895,7 @@ C This subroutine calculates multi-body contributions to hydrogen-bonding
 C Set lprn=.true. for debugging
       lprn=.false.
       eturn6=0.0d0
+      ecorr6=0.0d0
 #ifdef MPL
       n_corr=0
       n_corr1=0
@@ -5062,10 +5072,10 @@ cd                write(2,*)'wcorr6',wcorr6,' wturn6',wturn6
 cd                write(2,*)'ijkl',i,j,i+1,j1 
                 if (wcorr6.gt.0.0d0 .and. (j.ne.i+4 .or. j1.ne.i+3
      &               .or. wturn6.eq.0.0d0))then
-cd                  write (iout,*) '******ecorr6: i,j,i+1,j1',i,j,i+1,j1
-                  ecorr6=ecorr6+eello6(i,j,i+1,j1,jj,kk)
-cd                write (iout,*) 'ecorr',ecorr,' ecorr5=',ecorr5,
-cd     &            'ecorr6=',ecorr6
+c                  write (iout,*) '******ecorr6: i,j,i+1,j1',i,j,i+1,j1
+c                  ecorr6=ecorr6+eello6(i,j,i+1,j1,jj,kk)
+c                write (iout,*) 'ecorr',ecorr,' ecorr5=',ecorr5,
+c     &            'ecorr6=',ecorr6, wcorr6
 cd                write (iout,'(4e15.5)') sred_geom,
 cd     &          dabs(eello4(i,j,i+1,j1,jj,kk)),
 cd     &          dabs(eello5(i,j,i+1,j1,jj,kk)),
@@ -6388,7 +6398,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