Merge branch 'lipid' of mmka.chem.univ.gda.pl:unres into lipid
[unres.git] / source / cluster / wham / src / energy_p_new.F
index 622a1c6..959f6ba 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=iabs(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=iabs(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
@@ -2822,11 +2840,14 @@ 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. iabs(itype(iii)).eq.1 .and.
      &  iabs( 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)
@@ -2958,7 +2979,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,
       etheta=0.0D0
 c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
       do i=ithet_start,ithet_end
+         if (itype(i-1).eq.ntyp1) cycle
+         if (iabs(itype(i+1)).eq.20) iblock=2
+         if (iabs(itype(i+1)).ne.20) iblock=1
         dethetai=0.0d0
         dephii=0.0d0
         dephii1=0.0d0
         theti2=0.5d0*theta(i)
-CC Ta zmiana jest niewlasciwa
-        ityp2=ithetyp(iabs(itype(i-1)))
+        ityp2=ithetyp((itype(i-1)))
         do k=1,nntheterm
           coskt(k)=dcos(k*theti2)
           sinkt(k)=dsin(k*theti2)
@@ -3376,7 +3399,7 @@ CC Ta zmiana jest niewlasciwa
 #else
           phii=phi(i)
 #endif
-          ityp1=ithetyp(iabs(itype(i-2)))
+          ityp1=ithetyp((itype(i-2)))
           do k=1,nsingle
             cosph1(k)=dcos(k*phii)
             sinph1(k)=dsin(k*phii)
@@ -3397,7 +3420,7 @@ CC Ta zmiana jest niewlasciwa
 #else
           phii1=phi(i+1)
 #endif
-          ityp3=ithetyp(iabs(itype(i)))
+          ityp3=ithetyp((itype(i)))
           do k=1,nsingle
             cosph2(k)=dcos(k*phii1)
             sinph2(k)=dsin(k*phii1)
@@ -3413,7 +3436,7 @@ CC Ta zmiana jest niewlasciwa
 c        write (iout,*) "i",i," ityp1",itype(i-2),ityp1,
 c     &   " ityp2",itype(i-1),ityp2," ityp3",itype(i),ityp3
 c        call flush(iout)
-        ethetai=aa0thet(ityp1,ityp2,ityp3)
+        ethetai=aa0thet(ityp1,ityp2,ityp3,iblock)
         do k=1,ndouble
           do l=1,k-1
             ccl=cosph1(l)*cosph2(k-l)
@@ -3435,11 +3458,12 @@ c        call flush(iout)
         enddo
         endif
         do k=1,ntheterm
-          ethetai=ethetai+aathet(k,ityp1,ityp2,ityp3)*sinkt(k)
-          dethetai=dethetai+0.5d0*k*aathet(k,ityp1,ityp2,ityp3)
+          ethetai=ethetai+aathet(k,ityp1,ityp2,ityp3,iblock)*sinkt(k)
+          dethetai=dethetai+0.5d0*k*aathet(k,ityp1,ityp2,ityp3,iblock)
      &      *coskt(k)
           if (lprn)
-     &    write (iout,*) "k",k," aathet",aathet(k,ityp1,ityp2,ityp3),
+     &    write (iout,*) "k",k," aathet"
+     &       ,aathet(k,ityp1,ityp2,ityp3,iblock),
      &     " ethetai",ethetai
         enddo
         if (lprn) then
@@ -3458,24 +3482,24 @@ c        call flush(iout)
         endif
         do m=1,ntheterm2
           do k=1,nsingle
-            aux=bbthet(k,m,ityp1,ityp2,ityp3)*cosph1(k)
-     &         +ccthet(k,m,ityp1,ityp2,ityp3)*sinph1(k)
-     &         +ddthet(k,m,ityp1,ityp2,ityp3)*cosph2(k)
-     &         +eethet(k,m,ityp1,ityp2,ityp3)*sinph2(k)
+            aux=bbthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph1(k)
+     &         +ccthet(k,m,ityp1,ityp2,ityp3,iblock)*sinph1(k)
+     &         +ddthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph2(k)
+     &         +eethet(k,m,ityp1,ityp2,ityp3,iblock)*sinph2(k)
             ethetai=ethetai+sinkt(m)*aux
             dethetai=dethetai+0.5d0*m*aux*coskt(m)
             dephii=dephii+k*sinkt(m)*(
-     &          ccthet(k,m,ityp1,ityp2,ityp3)*cosph1(k)-
-     &          bbthet(k,m,ityp1,ityp2,ityp3)*sinph1(k))
+     &          ccthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph1(k)-
+     &          bbthet(k,m,ityp1,ityp2,ityp3,iblock)*sinph1(k))
             dephii1=dephii1+k*sinkt(m)*(
-     &          eethet(k,m,ityp1,ityp2,ityp3)*cosph2(k)-
-     &          ddthet(k,m,ityp1,ityp2,ityp3)*sinph2(k))
+     &          eethet(k,m,ityp1,ityp2,ityp3,iblock)*cosph2(k)-
+     &          ddthet(k,m,ityp1,ityp2,ityp3,iblock)*sinph2(k))
             if (lprn)
      &      write (iout,*) "m",m," k",k," bbthet",
-     &         bbthet(k,m,ityp1,ityp2,ityp3)," ccthet",
-     &         ccthet(k,m,ityp1,ityp2,ityp3)," ddthet",
-     &         ddthet(k,m,ityp1,ityp2,ityp3)," eethet",
-     &         eethet(k,m,ityp1,ityp2,ityp3)," ethetai",ethetai
+     &         bbthet(k,m,ityp1,ityp2,ityp3,iblock)," ccthet",
+     &         ccthet(k,m,ityp1,ityp2,ityp3,iblock)," ddthet",
+     &         ddthet(k,m,ityp1,ityp2,ityp3,iblock)," eethet",
+     &         eethet(k,m,ityp1,ityp2,ityp3,iblock)," ethetai",ethetai
           enddo
         enddo
         if (lprn)
@@ -3483,28 +3507,29 @@ c        call flush(iout)
         do m=1,ntheterm3
           do k=2,ndouble
             do l=1,k-1
-              aux=ffthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+
-     &            ffthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l)+
-     &            ggthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+
-     &            ggthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)
+              aux=ffthet(l,k,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(l,k)+
+     &            ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(k,l)+
+     &            ggthet(l,k,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(l,k)+
+     &            ggthet(k,l,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(k,l)
               ethetai=ethetai+sinkt(m)*aux
               dethetai=dethetai+0.5d0*m*coskt(m)*aux
               dephii=dephii+l*sinkt(m)*(
-     &           -ffthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)-
-     &            ffthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+
-     &            ggthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+
-     &            ggthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l))
+     &           -ffthet(l,k,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(l,k)-
+     &            ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(k,l)+
+     &            ggthet(l,k,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(l,k)+
+     &            ggthet(k,l,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(k,l))
               dephii1=dephii1+(k-l)*sinkt(m)*(
-     &           -ffthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+
-     &            ffthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+
-     &            ggthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)-
-     &            ggthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l))
+     &           -ffthet(l,k,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(l,k)+
+     &            ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(k,l)+
+     &            ggthet(l,k,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(l,k)-
+     &            ggthet(k,l,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(k,l))
               if (lprn) then
               write (iout,*) "m",m," k",k," l",l," ffthet",
-     &            ffthet(l,k,m,ityp1,ityp2,ityp3),
-     &            ffthet(k,l,m,ityp1,ityp2,ityp3)," ggthet",
-     &            ggthet(l,k,m,ityp1,ityp2,ityp3),
-     &            ggthet(k,l,m,ityp1,ityp2,ityp3)," ethetai",ethetai
+     &            ffthet(l,k,m,ityp1,ityp2,ityp3,iblock),
+     &            ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)," ggthet",
+     &            ggthet(l,k,m,ityp1,ityp2,ityp3,iblock),
+     &            ggthet(k,l,m,ityp1,ityp2,ityp3,iblock)," ethetai",
+     &            ethetai
               write (iout,*) cosph1ph2(l,k)*sinkt(m),
      &            cosph1ph2(k,l)*sinkt(m),
      &            sinph1ph2(l,k)*sinkt(m),sinph1ph2(k,l)*sinkt(m)
@@ -3867,7 +3892,7 @@ C     &   dc_norm(3,i+nres)
           y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
         enddo
         do j = 1,3
-          z_prime(j) = -uz(j,i-1)
+          z_prime(j) = -uz(j,i-1)*dsign(1.0d0,dfloat(itype(i)))
         enddo     
 c       write (2,*) "i",i
 c       write (2,*) "x_prime",(x_prime(j),j=1,3)
@@ -3890,7 +3915,7 @@ c
           xx = xx + x_prime(j)*dc_norm(j,i+nres)
           yy = yy + y_prime(j)*dc_norm(j,i+nres)
           zz = zz + z_prime(j)*dc_norm(j,i+nres)
-          zz = zz + dsign(1.0,itype(i))*z_prime(j)*dc_norm(j,i+nres)
+          zz = zz + z_prime(j)*dc_norm(j,i+nres)
         enddo
 
         xxtab(i)=xx
@@ -4080,8 +4105,10 @@ c     &   (dC_norm(j,i-1),j=1,3)," vbld_inv",vbld_inv(i+1),vbld_inv(i)
          dZZ_Ci1(k)=0.0d0
          dZZ_Ci(k)=0.0d0
          do j=1,3
-           dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1)*dC_norm(j,i+nres)
-           dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1)*dC_norm(j,i+nres)
+           dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1)
+     &      *dsign(1.0d0,dfloat(itype(i)))*dC_norm(j,i+nres)
+           dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1)
+     & *dsign(1.0d0,dfloat(itype(i)))*dC_norm(j,i+nres)
          enddo
           
          dXX_XYZ(k)=vbld_inv(i+nres)*(x_prime(k)-xx*dC_norm(k,i+nres))
@@ -4904,6 +4931,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
@@ -5080,10 +5108,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)),
@@ -6406,7 +6434,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