removed uninitialized ghalf from equation
[unres.git] / source / unres / src_MD-M / energy_p_new_barrier.F
index 5ad53f1..700fb07 100644 (file)
@@ -1420,6 +1420,9 @@ C
       include 'COMMON.CONTROL'
       include 'COMMON.SBRIDGE'
       logical lprn
+
+c      write(iout,*) "Jestem w egb(evdw)"
+
       evdw=0.0D0
 ccccc      energy_dec=.false.
 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
         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
+
+c              write(iout,*) "PRZED ZWYKLE", evdwij
               call dyn_ssbond_ene(i,j,evdwij)
+c              write(iout,*) "PO ZWYKLE", evdwij
+
               evdw=evdw+evdwij
               if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') 
      &                        'evdw',i,j,evdwij,' ss'
+C triple bond artifac removal
+             do k=j+1,iend(i,iint) 
+C search over all next residues
+              if (dyn_ss_mask(k)) then
+C check if they are cysteins
+C              write(iout,*) 'k=',k
+
+c              write(iout,*) "PRZED TRI", evdwij
+               evdwij_przed_tri=evdwij
+              call triple_ssbond_ene(i,j,k,evdwij)
+c               if(evdwij_przed_tri.ne.evdwij) then
+c                 write (iout,*) "TRI:", evdwij, evdwij_przed_tri
+c               endif
+
+c              write(iout,*) "PO TRI", evdwij
+C call the energy function that removes the artifical triple disulfide
+C bond the soubroutine is located in ssMD.F
+              evdw=evdw+evdwij             
+              if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
+     &                        'evdw',i,j,evdwij,'tss'
+              endif!dyn_ss_mask(k)
+             enddo! k
             ELSE
             ind=ind+1
             itypj=iabs(itype(j))
@@ -4085,44 +4114,44 @@ C 15/02/13 CC dynamic SSbond - additional check
      & iabs(itype(jjj)).eq.1) then
           call ssbond_ene(iii,jjj,eij)
           ehpb=ehpb+2*eij
-cd          write (iout,*) "eij",eij
          endif
+cd          write (iout,*) "eij",eij
         else
 C Calculate the distance between the two points and its difference from the
 C target distance.
-        dd=dist(ii,jj)
-        rdis=dd-dhpb(i)
+          dd=dist(ii,jj)
+            rdis=dd-dhpb(i)
 C Get the force constant corresponding to this distance.
-        waga=forcon(i)
+            waga=forcon(i)
 C Calculate the contribution to energy.
-        ehpb=ehpb+waga*rdis*rdis
+            ehpb=ehpb+waga*rdis*rdis
 C
 C Evaluate gradient.
 C
-        fac=waga*rdis/dd
+            fac=waga*rdis/dd
 cd      print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd,
 cd   &   ' waga=',waga,' fac=',fac
-        do j=1,3
-          ggg(j)=fac*(c(j,jj)-c(j,ii))
-        enddo
+            do j=1,3
+              ggg(j)=fac*(c(j,jj)-c(j,ii))
+            enddo
 cd      print '(i3,3(1pe14.5))',i,(ggg(j),j=1,3)
 C If this is a SC-SC distance, we need to calculate the contributions to the
 C Cartesian gradient in the SC vectors (ghpbx).
-        if (iii.lt.ii) then
+          if (iii.lt.ii) then
           do j=1,3
             ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
             ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
           enddo
-        endif
+          endif
 cgrad        do j=iii,jjj-1
 cgrad          do k=1,3
 cgrad            ghpbc(k,j)=ghpbc(k,j)+ggg(k)
 cgrad          enddo
 cgrad        enddo
-        do k=1,3
-          ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
-          ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
-        enddo
+          do k=1,3
+            ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
+            ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
+          enddo
         endif
       enddo
       ehpb=0.5D0*ehpb
@@ -4574,7 +4603,9 @@ C
       logical lprn /.false./, lprn1 /.false./
       etheta=0.0D0
       do i=ithet_start,ithet_end
-        if (itype(i-1).eq.ntyp1) cycle
+        if ((itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or.
+     &(itype(i).eq.ntyp1)) cycle
+C        print *,i,theta(i)
         if (iabs(itype(i+1)).eq.20) iblock=2
         if (iabs(itype(i+1)).ne.20) iblock=1
         dethetai=0.0d0
@@ -4586,7 +4617,9 @@ C
           coskt(k)=dcos(k*theti2)
           sinkt(k)=dsin(k*theti2)
         enddo
-        if (i.gt.3 .and. itype(i-2).ne.ntyp1) then
+C        print *,ethetai
+
+        if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
 #ifdef OSF
           phii=phi(i)
           if (phii.ne.phii) phii=150.0
@@ -4601,13 +4634,13 @@ C propagation of chirality for glycine type
           enddo
         else
           phii=0.0d0
-          ityp1=nthetyp+1
           do k=1,nsingle
+          ityp1=ithetyp((itype(i-2)))
             cosph1(k)=0.0d0
             sinph1(k)=0.0d0
           enddo 
         endif
-        if (i.lt.nres .and. itype(i).ne.ntyp1) 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
@@ -4622,7 +4655,7 @@ C propagation of chirality for glycine type
           enddo
         else
           phii1=0.0d0
-          ityp3=nthetyp+1
+          ityp3=ithetyp((itype(i)))
           do k=1,nsingle
             cosph2(k)=0.0d0
             sinph2(k)=0.0d0
@@ -4672,6 +4705,7 @@ C propagation of chirality for glycine type
         enddo
         write(iout,*) "ethetai",ethetai
         endif
+C       print *,ethetai
         do m=1,ntheterm2
           do k=1,nsingle
             aux=bbthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph1(k)
@@ -4692,10 +4726,16 @@ C propagation of chirality for glycine type
      &         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
+C        print *,"tu",cosph1(k),sinph1(k),cosph2(k),sinph2(k)
           enddo
         enddo
+C        print *,"cosph1", (cosph1(k), k=1,nsingle)
+C        print *,"cosph2", (cosph2(k), k=1,nsingle)
+C        print *,"sinph1", (sinph1(k), k=1,nsingle)
+C        print *,"sinph2", (sinph2(k), k=1,nsingle)
         if (lprn)
      &  write(iout,*) "ethetai",ethetai
+C        print *,"tu",cosph1(k),sinph1(k),cosph2(k),sinph2(k)
         do m=1,ntheterm3
           do k=2,ndouble
             do l=1,k-1
@@ -4731,6 +4771,7 @@ C propagation of chirality for glycine type
         enddo
 10      continue
 c        lprn1=.true.
+C        print *,ethetai
         if (lprn1) 
      &    write (iout,'(i2,3f8.1,9h ethetai ,f10.5)') 
      &   i,theta(i)*rad2deg,phii*rad2deg,
@@ -4739,7 +4780,7 @@ 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
+        gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*dethetai
       enddo
       return
       end
@@ -7733,9 +7774,9 @@ cd        ghalf=0.0d0
 cold        ghalf=0.5d0*eel5*eij*gacont_hbr(ll,kk,k)
 cgrad        ghalf=0.5d0*ggg2(ll)
 cd        ghalf=0.0d0
-        gradcorr5(ll,k)=gradcorr5(ll,k)+ghalf+ekont*derx(ll,2,2)
+        gradcorr5(ll,k)=gradcorr5(ll,k)+ekont*derx(ll,2,2)
         gradcorr5(ll,k+1)=gradcorr5(ll,k+1)+ekont*derx(ll,3,2)
-        gradcorr5(ll,l)=gradcorr5(ll,l)+ghalf+ekont*derx(ll,4,2)
+        gradcorr5(ll,l)=gradcorr5(ll,l)+ekont*derx(ll,4,2)
         gradcorr5(ll,l1)=gradcorr5(ll,l1)+ekont*derx(ll,5,2)
         gradcorr5_long(ll,l)=gradcorr5_long(ll,l)+gradcorr5kl
         gradcorr5_long(ll,k)=gradcorr5_long(ll,k)-gradcorr5kl
@@ -8040,12 +8081,12 @@ C                                                                              C
 C          o             o                                                     C
 C     \   /l\           /j\   /                                                C
 C      \ /   \         /   \ /                                                 C
-C       o| o |         | o |o                                                  C
+C       o| o |         | o |o                                                  C                
 C     \ j|/k\|      \  |/k\|l                                                  C
 C      \ /   \       \ /   \                                                   C
 C       o             o                                                        C
-C       i             i                                                        C
-C                                                                              C
+C       i             i                                                        C 
+C                                                                              C           
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 cd      write (2,*) 'eello6_graph2: i,',i,' j',j,' k',k,' l',l
 C AL 7/4/01 s1 would occur in the sixth-order moment, 
@@ -8216,10 +8257,10 @@ c----------------------------------------------------------------------------
       double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2)
       logical swap
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-C                                                                              C
+C                                                                              C 
 C      Parallel       Antiparallel                                             C
 C                                                                              C
-C          o             o                                                     C
+C          o             o                                                     C 
 C         /l\   /   \   /j\                                                    C 
 C        /   \ /     \ /   \                                                   C
 C       /| o |o       o| o |\                                                  C
@@ -8333,7 +8374,7 @@ c----------------------------------------------------------------------------
      & auxvec1(2),auxmat1(2,2)
       logical swap
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-C                                                                              C
+C                                                                              C                       
 C      Parallel       Antiparallel                                             C
 C                                                                              C
 C          o             o                                                     C
@@ -8341,10 +8382,10 @@ C         /l\   /   \   /j\                                                    C
 C        /   \ /     \ /   \                                                   C
 C       /| o |o       o| o |\                                                  C
 C     \ j|/k\|      \  |/k\|l                                                  C
-C      \ /   \       \ /   \                                                   C
+C      \ /   \       \ /   \                                                   C 
 C       o     \       o     \                                                  C
 C       i             i                                                        C
-C                                                                              C
+C                                                                              C 
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 C
 C 4/7/01 AL Component s1 was removed, because it pertains to the respective