The wrong one gradient compontent...

 cleaning in progress by A.Liwo other gradient OK

diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F
index 8f89ae9..e56e104 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -2404,13 +2404,13 @@ c        if (i.gt. iatel_s+2 .and. i.lt.iatel_e+5) then
         if (i.gt. nnt+2 .and. i.lt.nct+2) then
           iti = itortyp(itype(i-2))
         else
-          iti=ntortyp+1
+          iti=ntortyp
         endif
 c        if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
         if (i.gt. nnt+1 .and. i.lt.nct+1) then
           iti1 = itortyp(itype(i-1))
         else
-          iti1=ntortyp+1
+          iti1=ntortyp
         endif
 cd        write (iout,*) '*******i',i,' iti1',iti
 cd        write (iout,*) 'b1',b1(:,iti)
@@ -2451,10 +2451,10 @@ c        if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
           if (itype(i-1).le.ntyp) then
             iti1 = itortyp(itype(i-1))
           else
-            iti1=ntortyp+1
+            iti1=ntortyp
           endif
         else
-          iti1=ntortyp+1
+          iti1=ntortyp
         endif
         do k=1,2
           mu(k,i-2)=Ub2(k,i-2)+b1(k,iti1)
@@ -2869,7 +2869,11 @@ C
 C 14/01/2014 TURN3,TUNR4 does no go under periodic boundry condition
       do i=iturn3_start,iturn3_end
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
-     &  .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1) cycle
+     &  .or. itype(i+2).eq.ntyp1
+     &  .or. itype(i+3).eq.ntyp1
+c     &  .or. itype(i-1).eq.ntyp1
+     &  .or. itype(i+4).eq.ntyp1
+     &  ) cycle
         dxi=dc(1,i)
         dyi=dc(2,i)
         dzi=dc(3,i)
@@ -2912,7 +2916,9 @@ C Condition for being inside the proper box
       do i=iturn4_start,iturn4_end
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
      &    .or. itype(i+3).eq.ntyp1
-     &    .or. itype(i+4).eq.ntyp1) cycle
+     &    .or. itype(i+4).eq.ntyp1
+     &    .or. itype(i+5).eq.ntyp1
+     &                             ) cycle
         dxi=dc(1,i)
         dyi=dc(2,i)
         dzi=dc(3,i)
@@ -2962,7 +2968,9 @@ c
 c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
 c
       do i=iatel_s,iatel_e
-        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
+     &  .or. itype(i+2).eq.ntyp1
+     &                ) cycle
         dxi=dc(1,i)
         dyi=dc(2,i)
         dzi=dc(3,i)
@@ -3002,7 +3010,9 @@ c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
         num_conti=num_cont_hb(i)
         do j=ielstart(i),ielend(i)
 c          write (iout,*) i,j,itype(i),itype(j)
-          if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1) cycle
+          if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1
+     & .or.itype(j+2).eq.ntyp1
+     &) cycle
           call eelecij(i,j,ees,evdw1,eel_loc)
         enddo ! j
         num_cont_hb(i)=num_conti
@@ -3186,9 +3196,15 @@ cgrad            do l=1,3
 cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
 cgrad            enddo
 cgrad          enddo
-          ggg(1)=facvdw*xj
-          ggg(2)=facvdw*yj
-          ggg(3)=facvdw*zj
+          if (sss.gt.0.0) then
+          ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj
+          ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj
+          ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj
+          else
+          ggg(1)=0.0
+          ggg(2)=0.0
+          ggg(3)=0.0
+          endif
 c          do k=1,3
 c            ghalf=0.5D0*ggg(k)
 c            gvdwpp(k,i)=gvdwpp(k,i)+ghalf
@@ -3212,7 +3228,7 @@ C MARYSIA
           facvdw=(ev1+evdwij)*sss
           facel=(el1+eesij)
           fac1=fac
-          fac=-3*rrmij*(facvdw+facvdw+facel)+sssgrad*rmij*evdwij
+          fac=-3*rrmij*(facvdw+facvdw+facel)
           erij(1)=xj*rmij
           erij(2)=yj*rmij
           erij(3)=zj*rmij
@@ -3241,9 +3257,9 @@ cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
 cgrad            enddo
 cgrad          enddo
 c 9/28/08 AL Gradient compotents will be summed only at the end
-          ggg(1)=facvdw*xj*sss
-          ggg(2)=facvdw*yj*sss
-          ggg(3)=facvdw*zj*sss
+          ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj
+          ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj
+          ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj
           do k=1,3
             gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
             gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
@@ -3479,7 +3495,8 @@ cgrad            endif
 C Contribution to the local-electrostatic energy coming from the i-j pair
           eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
      &     +a33*muij(4)
-cd          write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
+c          write (iout,*) 'i',i,' j',j,itype(i),itype(j),
+c     &                     ' eel_loc_ij',eel_loc_ij
 
           if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
      &            'eelloc',i,j,eel_loc_ij
@@ -4652,9 +4669,8 @@ c      time12=1.0d0
       etheta=0.0D0
 c     write (*,'(a,i2)') 'EBEND ICG=',icg
       do i=ithet_start,ithet_end
-        print *,i,itype(i-1),itype(i),itype(i-2)
-        if (itype(i-1).eq.ntyp1) cycle
-        print *,'wchodze',itype(i-1)
+        if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
+     &  .or.itype(i).eq.ntyp1) cycle
 C Zero the energy function and its derivative at 0 or pi.
         call splinthet(theta(i),0.5d0*delta,ss,ssd)
         it=itype(i-1)
@@ -4671,7 +4687,7 @@ C Zero the energy function and its derivative at 0 or pi.
           ichir22=isign(1,itype(i))
          endif
 
-        if (i.gt.3 .and. itype(i-2).ne.ntyp1) then
+        if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
 #ifdef OSF
          phii=phi(i)
           if (phii.ne.phii) phii=150.0
@@ -4684,7 +4700,7 @@ C Zero the energy function and its derivative at 0 or pi.
           y(1)=0.0D0
           y(2)=0.0D0
         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
@@ -4692,8 +4708,8 @@ C Zero the energy function and its derivative at 0 or pi.
           z(1)=cos(phii1)
 #else
           phii1=phi(i+1)
-          z(1)=dcos(phii1)
 #endif
+          z(1)=dcos(phii1)
           z(2)=dsin(phii1)
         else
           z(1)=0.0D0
@@ -4897,7 +4913,10 @@ C
       etheta=0.0D0
       do i=ithet_start,ithet_end
 c        print *,i,itype(i-1),itype(i),itype(i-2)
-        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 In current verion the ALL DUMMY ATOM POTENTIALS ARE OFF
+
         if (iabs(itype(i+1)).eq.20) iblock=2
         if (iabs(itype(i+1)).ne.20) iblock=1
         dethetai=0.0d0
@@ -4909,7 +4928,7 @@ c        print *,i,itype(i-1),itype(i),itype(i-2)
           coskt(k)=dcos(k*theti2)
           sinkt(k)=dsin(k*theti2)
         enddo
-        if (i.gt.3 .and. itype(i-2).ne.ntyp1) then
+        if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
 #ifdef OSF
           phii=phi(i)
           if (phii.ne.phii) phii=150.0
@@ -4930,7 +4949,7 @@ C propagation of chirality for glycine type
             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
@@ -5828,9 +5847,9 @@ c      lprn=.true.
       do i=iphi_start,iphi_end
       etors_ii=0.0D0
         if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1
-     &      .or. itype(i).eq.ntyp1) cycle
-       itori=itortyp(itype(i-2))
-       itori1=itortyp(itype(i-1))
+     &      .or. itype(i).eq.ntyp1 .or. itype(i-3).eq.ntyp1) cycle
+        itori=itortyp(itype(i-2))
+        itori1=itortyp(itype(i-1))
         phii=phi(i)
         gloci=0.0D0
 C Proline-Proline pair is a special case...
@@ -5928,8 +5947,9 @@ C ANY TWO ARE DUMMY ATOMS in row CYCLE
 c        if (((itype(i-3).eq.ntyp1).and.(itype(i-2).eq.ntyp1)).or.
 c     &      ((itype(i-2).eq.ntyp1).and.(itype(i-1).eq.ntyp1))  .or.
 c     &      ((itype(i-1).eq.ntyp1).and.(itype(i).eq.ntyp1))) cycle
-         if ((itype(i-3).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or.
-     &  (itype(i-1).eq.ntyp1).or.(itype(i).eq.ntyp1)) cycle
+        if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1
+     &      .or. itype(i).eq.ntyp1 .or. itype(i-3).eq.ntyp1) cycle
+C In current verion the ALL DUMMY ATOM POTENTIALS ARE OFF
 C For introducing the NH3+ and COO- group please check the etor_d for reference
 C and guidance
         etors_ii=0.0D0
@@ -6033,10 +6053,14 @@ c     lprn=.true.
 c      write(iout,*) "a tu??"
       do i=iphid_start,iphid_end
 C ANY TWO ARE DUMMY ATOMS in row CYCLE
-        if (((itype(i-3).eq.ntyp1).and.(itype(i-2).eq.ntyp1)).or.
-     &      ((itype(i-2).eq.ntyp1).and.(itype(i-1).eq.ntyp1)).or.
-     &      ((itype(i-1).eq.ntyp1).and.(itype(i).eq.ntyp1))  .or.
-     &      ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1))) cycle
+C        if (((itype(i-3).eq.ntyp1).and.(itype(i-2).eq.ntyp1)).or.
+C     &      ((itype(i-2).eq.ntyp1).and.(itype(i-1).eq.ntyp1)).or.
+C     &      ((itype(i-1).eq.ntyp1).and.(itype(i).eq.ntyp1))  .or.
+C     &      ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1))) cycle
+         if ((itype(i-2).eq.ntyp1).or.itype(i-3).eq.ntyp1.or.
+     &  (itype(i-1).eq.ntyp1).or.(itype(i).eq.ntyp1).or.
+     &  (itype(i+1).eq.ntyp1)) cycle
+C In current verion the ALL DUMMY ATOM POTENTIALS ARE OFF
         itori=itortyp(itype(i-2))
         itori1=itortyp(itype(i-1))
         itori2=itortyp(itype(i))
@@ -6046,6 +6070,7 @@ C ANY TWO ARE DUMMY ATOMS in row CYCLE
         gloci2=0.0D0
         iblock=1
         if (iabs(itype(i+1)).eq.20) iblock=2
+C Iblock=2 Proline type
 C ADASKO: WHEN PARAMETERS FOR THIS TYPE OF BLOCKING GROUP IS READY UNCOMMENT
 C CHECK WEATHER THERE IS NECCESITY FOR iblock=3 for COO-
 C        if (itype(i+1).eq.ntyp1) iblock=3
@@ -7174,7 +7199,7 @@ C---------------------------------------------------------------------------
       if (j.lt.nres-1) then
         itj1 = itortyp(itype(j+1))
       else
-        itj1=ntortyp+1
+        itj1=ntortyp
       endif
       do iii=1,2
         dipi(iii,1)=Ub2(iii,i)
@@ -7264,14 +7289,14 @@ C parallel orientation of the two CA-CA-CA frames.
         if (i.gt.1) then
           iti=itortyp(itype(i))
         else
-          iti=ntortyp+1
+          iti=ntortyp
         endif
         itk1=itortyp(itype(k+1))
         itj=itortyp(itype(j))
         if (l.lt.nres-1) then
           itl1=itortyp(itype(l+1))
         else
-          itl1=ntortyp+1
+          itl1=ntortyp
         endif
 C A1 kernel(j+1) A2T
 cd        do iii=1,2
@@ -7417,7 +7442,7 @@ C Antiparallel orientation of the two CA-CA-CA frames.
         if (i.gt.1) then
           iti=itortyp(itype(i))
         else
-          iti=ntortyp+1
+          iti=ntortyp
         endif
         itk1=itortyp(itype(k+1))
         itl=itortyp(itype(l))
@@ -7425,7 +7450,7 @@ C Antiparallel orientation of the two CA-CA-CA frames.
         if (j.lt.nres-1) then
           itj1=itortyp(itype(j+1))
         else 
-          itj1=ntortyp+1
+          itj1=ntortyp
         endif
 C A2 kernel(j-1)T A1T
         call kernel(aa1(1,1),aa2t(1,1),a_chuj_der(1,1,1,1,jj,i),
@@ -8579,14 +8604,14 @@ C           energy moment and not to the cluster cumulant.
       if (j.lt.nres-1) then
         itj1=itortyp(itype(j+1))
       else
-        itj1=ntortyp+1
+        itj1=ntortyp
       endif
       itk=itortyp(itype(k))
       itk1=itortyp(itype(k+1))
       if (l.lt.nres-1) then
         itl1=itortyp(itype(l+1))
       else
-        itl1=ntortyp+1
+        itl1=ntortyp
       endif
 #ifdef MOMENT
       s1=dip(4,jj,i)*dip(4,kk,k)
@@ -8698,19 +8723,19 @@ cd      write (2,*) 'eello_graph4: wturn6',wturn6
       if (j.lt.nres-1) then
         itj1=itortyp(itype(j+1))
       else
-        itj1=ntortyp+1
+        itj1=ntortyp
       endif
       itk=itortyp(itype(k))
       if (k.lt.nres-1) then
         itk1=itortyp(itype(k+1))
       else
-        itk1=ntortyp+1
+        itk1=ntortyp
       endif
       itl=itortyp(itype(l))
       if (l.lt.nres-1) then
         itl1=itortyp(itype(l+1))
       else
-        itl1=ntortyp+1
+        itl1=ntortyp
       endif
 cd      write (2,*) 'eello6_graph4:','i',i,' j',j,' k',k,' l',l
 cd      write (2,*) 'iti',iti,' itj',itj,' itj1',itj1,' itk',itk,

diff --git a/source/unres/src_MD-M/initialize_p.F b/source/unres/src_MD-M/initialize_p.F
index 70e35ab..1a5e34a 100644 (file)
--- a/source/unres/src_MD-M/initialize_p.F
+++ b/source/unres/src_MD-M/initialize_p.F
@@ -245,6 +245,22 @@ C Initialize the bridge arrays
        ihpb(i)=0
        jhpb(i)=0
       enddo
+C Initialize correlation arrays
+      do i=-maxtor,maxtor
+       do k=1,2
+        b1(k,i)=0.0
+        b2(k,i)=0.0
+        b1tilde(k,i)=0.0
+c        b2tilde(k,i)=0.0
+        do j=1,2
+        CC(j,k,i)=0.0
+        Ctilde(j,k,i)=0.0
+        DD(j,k,i)=0.0
+        Dtilde(j,k,i)=0.0
+        EE(j,k,i)=0.0
+        enddo
+       enddo
+      enddo
 C
 C Initialize timing.
 C

diff --git a/source/unres/src_MD-M/intcartderiv.F b/source/unres/src_MD-M/intcartderiv.F
index d715c19..12d7dee 100644 (file)
--- a/source/unres/src_MD-M/intcartderiv.F
+++ b/source/unres/src_MD-M/intcartderiv.F
@@ -99,7 +99,8 @@ c conventional formulas around 0 and 180.
 #else
       do i=4,nres      
 #endif
-c        if (itype(i-1).eq.21 .or. itype(i-2).eq.21 ) cycle
+c        if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1
+c     &      .or. itype(i).eq.ntyp1 .or. itype(i-3).eq.ntyp1) cycle
 c the conventional case
         sint=dsin(theta(i))
        sint1=dsin(theta(i-1))
@@ -124,8 +125,8 @@ c    Obtaining the gamma derivatives from sine derivative
             ctgt=cost/sint
             ctgt1=cost1/sint1
             cosg_inv=1.0d0/cosg
-            if (itype(i-1).ne.ntyp1 .and. itype(i-2).ne.ntyp1) then
-           dsinphi(j,1,i)=-sing*ctgt1*dtheta(j,1,i-1)
+c            if (itype(i-1).ne.ntyp1 .and. itype(i-2).ne.ntyp1) then
+      dsinphi(j,1,i)=-sing*ctgt1*dtheta(j,1,i-1)
      &        -(fac0*vp1(j)+sing*dc_norm(j,i-3))*vbld_inv(i-2)
             dphi(j,1,i)=cosg_inv*dsinphi(j,1,i)
             dsinphi(j,2,i)=
@@ -136,13 +137,13 @@ c    Obtaining the gamma derivatives from sine derivative
      &        +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i)
 c     &        +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
             dphi(j,3,i)=cosg_inv*dsinphi(j,3,i)
-            endif
+c            endif
 c Bug fixed 3/24/05 (AL)
         enddo                                              
 c   Obtaining the gamma derivatives from cosine derivative
         else
            do j=1,3
-           if (itype(i-1).ne.ntyp1 .and. itype(i-2).ne.ntyp1) then
+c           if (itype(i-1).ne.ntyp1 .and. itype(i-2).ne.ntyp1) then
            dcosphi(j,1,i)=fac1*dcostheta(j,1,i-1)+fac3*
      &    dcostheta(j,1,i-1)-fac0*(dc_norm(j,i-1)-scalp*
      &     dc_norm(j,i-3))/vbld(i-2)
@@ -155,7 +156,7 @@ c   Obtaining the gamma derivatives from cosine derivative
      &    dcostheta(j,2,i)-fac0*(dc_norm(j,i-3)-scalp*
      &     dc_norm(j,i-1))/vbld(i)
            dphi(j,3,i)=-1/sing*dcosphi(j,3,i)       
-           endif
+c           endif
          enddo
         endif                                                                                           
       enddo

diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F
index 22236bc..bd2165b 100644 (file)
--- a/source/unres/src_MD-M/parmread.F
+++ b/source/unres/src_MD-M/parmread.F
@@ -591,7 +591,6 @@ C
       read (itorp,*,end=113,err=113) ntortyp,nterm_old
       if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
-      itortyp(ntyp1)=0
       do i=1,ntortyp
        do j=1,ntortyp
          read (itorp,'(a)')