Dzialajacy gradient dla reszt 13 i 7

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 c13d341..2682840 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -2252,7 +2252,7 @@ C
 C Compute the virtual-bond-torsional-angle dependent quantities needed
 C to calculate the el-loc multibody terms of various order.
 C
-      write(iout,*) 'nphi=',nphi,nres
+c      write(iout,*) 'nphi=',nphi,nres
 #ifdef PARMAT
       do i=ivec_start+2,ivec_end+2
 #else
@@ -2270,23 +2270,23 @@ c        if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
         else
           iti1=ntortyp+1
         endif
-        write(iout,*),i
-        b1(1,i-2)=bnew1(1,1,iti)*sin(theta(i-1)/2.0)
-     &           +bnew1(2,1,iti)*sin(theta(i-1))
-     &           +bnew1(3,1,iti)*cos(theta(i-1)/2.0)
-        gtb1(1,i-2)=bnew1(1,1,iti)*cos(theta(i-1)/2.0)/2.0
-     &             +bnew1(2,1,iti)*cos(theta(i-1))
-     &             -bnew1(3,1,iti)*sin(theta(i-1)/2.0)/2.0
+c        write(iout,*),i
+        b1(1,i-2)=bnew1(1,1,iti)*dsin(theta(i-1)/2.0)
+     &           +bnew1(2,1,iti)*dsin(theta(i-1))
+     &           +bnew1(3,1,iti)*dcos(theta(i-1)/2.0)
+        gtb1(1,i-2)=bnew1(1,1,iti)*dcos(theta(i-1)/2.0d0)/2.0d0
+     &             +bnew1(2,1,iti)*dcos(theta(i-1))
+     &             -bnew1(3,1,iti)*dsin(theta(i-1)/2.0d0)/2.0d0
 c     &           +bnew1(3,1,iti)*sin(alpha(i))*cos(beta(i))
 c     &*(cos(theta(i)/2.0)
-        b2(1,i-2)=bnew2(1,1,iti)*sin(theta(i-1)/2.0)
-     &           +bnew2(2,1,iti)*sin(theta(i-1))
-     &           +bnew2(3,1,iti)*cos(theta(i-1)/2.0)
+        b2(1,i-2)=bnew2(1,1,iti)*dsin(theta(i-1)/2.0)
+     &           +bnew2(2,1,iti)*dsin(theta(i-1))
+     &           +bnew2(3,1,iti)*dcos(theta(i-1)/2.0)
 c     &           +bnew2(3,1,iti)*sin(alpha(i))*cos(beta(i))
 c     &*(cos(theta(i)/2.0)
-        gtb2(1,i-2)=bnew2(1,1,iti)*cos(theta(i-1)/2.0)/2.0
-     &             +bnew2(2,1,iti)*cos(theta(i-1))
-     &             -bnew2(3,1,iti)*sin(theta(i-1)/2.0)/2.0
+        gtb2(1,i-2)=bnew2(1,1,iti)*dcos(theta(i-1)/2.0d0)/2.0d0
+     &             +bnew2(2,1,iti)*dcos(theta(i-1))
+     &             -bnew2(3,1,iti)*dsin(theta(i-1)/2.0d0)/2.0d0
 c        if (ggb1(1,i).eq.0.0d0) then
 c        write(iout,*) 'i=',i,ggb1(1,i),
 c     &bnew1(1,1,iti)*cos(theta(i)/2.0)/2.0,
@@ -2307,8 +2307,9 @@ c        b2(2,iti)=bnew2(1,2,iti)*sin(alpha(i))*sin(beta(i))
        b1tilde(2,i-2)=-b1(2,i-2)
        b2tilde(1,i-2)=b2(1,i-2)
        b2tilde(2,i-2)=-b2(2,i-2)
-       write (iout,*) 'i=',i-2,gtb1(2,i-2),gtb1(1,i-2)
-       write (iout,*) 'theta=', theta(i-1)
+c       write (iout,*) 'i=',i-2,gtb1(2,i-2),gtb1(1,i-2)
+c       write(iout,*)  'b1=',b1(1,i-2)
+c       write (iout,*) 'theta=', theta(i-1)
        enddo
 #ifdef PARMAT
       do i=ivec_start+2,ivec_end+2
@@ -2447,7 +2448,7 @@ c        if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
         do k=1,2
           mu(k,i-2)=Ub2(k,i-2)+b1(k,i-1)
         enddo
-cd        write (iout,*) 'mu ',mu(:,i-2)
+c        write (iout,*) 'mu ',mu(:,i-2),i-2
 cd        write (iout,*) 'mu1',mu1(:,i-2)
 cd        write (iout,*) 'mu2',mu2(:,i-2)
         if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.wcorr6.gt.0.0d0)
@@ -2895,7 +2896,7 @@ c
 c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
 c
 c      do i=iatel_s,iatel_e
-       do i=4,5
+       do i=7,7
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
         dxi=dc(1,i)
         dyi=dc(2,i)
@@ -2909,8 +2910,8 @@ c      do i=iatel_s,iatel_e
 c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
         num_conti=num_cont_hb(i)
 c        do j=ielstart(i),ielend(i)
-         do j=8,9
-          write (iout,*) 'tu wchodze',i,j,itype(i),itype(j)
+         do j=13,13
+c          write (iout,*) 'tu wchodze',i,j,itype(i),itype(j)
           if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1) cycle
           call eelecij(i,j,ees,evdw1,eel_loc)
         enddo ! j
@@ -3181,13 +3182,14 @@ C
             do l=1,2
               kkk=kkk+1
               muij(kkk)=mu(k,i)*mu(l,j)
+c              write(iout,*) 'mumu=', mu(k,i),mu(l,j),i,j,k,l
 #ifdef NEWCORR
-             gmuij1(kkk)=gtb1(k,i)*mu(l,j)
-             write(iout,*) 'kkk=', gtb1(k,i)*mu(l,j),gtb1(k,i),k,i
-             gmuij2(kkk)=gUb2(k,i-1)*mu(l,j)
-             gmuji1(kkk)=mu(k,i)*gtb1(l,j)
-             write(iout,*) 'kkk=', gtb1(k,i)*mu(l,j),gtb1(l,j),l,j
-             gmuji2(kkk)=mu(k,i)*gUb2(l,j-1)
+             gmuij1(kkk)=gtb1(k,i+1)*mu(l,j)
+c             write(iout,*) 'k=',k,i,gtb1(k,i+1),gtb1(k,i+1)*mu(l,j)
+             gmuij2(kkk)=gUb2(k,i)*mu(l,j)
+             gmuji1(kkk)=mu(k,i)*gtb1(l,j+1)
+c             write(iout,*) 'l=',l,j,gtb1(l,j+1),gtb1(l,j+1)*mu(k,i)
+             gmuji2(kkk)=mu(k,i)*gUb2(l,j)
 #endif
             enddo
           enddo  
@@ -3354,37 +3356,47 @@ 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)
+c          write(iout,*) 'muije=',muij(1),muij(2),muij(3),muij(4)
 C Calculate patrial derivative for theta angle
 #ifdef NEWCORR
          geel_loc_ij=a22*gmuij1(1)
      &     +a23*gmuij1(2)
      &     +a32*gmuij1(3)
      &     +a33*gmuij1(4)         
-         write(iout,*) "derivative over thatai"
-         write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
-     &   a33*gmuij1(4) 
+c         write(iout,*) "derivative over thatai"
+c         write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
+c     &   a33*gmuij1(4) 
          gloc(nphi+i,icg)=gloc(nphi+i,icg)+
      &      geel_loc_ij*wel_loc
-         write(iout,*) "derivative over thatai-1" 
-         write(iout,*) a22*gmuij2(1), a23*gmuij2(2) ,a32*gmuij2(3),
-     &   a33*gmuij2(4)
-         geel_loc_ij=a22*gmuij2(1)+a23*gmuij2(2)+a32*gmuij2(3)
+c         write(iout,*) "derivative over thatai-1" 
+c         write(iout,*) a22*gmuij2(1), a23*gmuij2(2) ,a32*gmuij2(3),
+c     &   a33*gmuij2(4)
+         geel_loc_ij=
+     &     a22*gmuij2(1)
+     &     +a23*gmuij2(2)
+     &     +a32*gmuij2(3)
      &     +a33*gmuij2(4)
          gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
      &      geel_loc_ij*wel_loc
-         geel_loc_ji=a22*gmuji1(1)+a23*gmuji1(2)+a32*gmuji1(3)
+c  Derivative over j residue
+         geel_loc_ji=a22*gmuji1(1)
+     &     +a23*gmuji1(2)
+     &     +a32*gmuji1(3)
      &     +a33*gmuji1(4)
-         write(iout,*) "derivative over thataj" 
-         write(iout,*) a22*gmuji1(1), a23*gmuji1(2) ,a32*gmuji1(3),
-     &   a33*gmuji1(4)
+c         write(iout,*) "derivative over thataj" 
+c         write(iout,*) a22*gmuji1(1), a23*gmuji1(2) ,a32*gmuji1(3),
+c     &   a33*gmuji1(4)
 
-         gloc(nphi+j,icg)=gloc(nphi+j,icg)+
+        gloc(nphi+j,icg)=gloc(nphi+j,icg)+
      &      geel_loc_ji*wel_loc
-         geel_loc_ji=a22*gmuji2(1)+a23*gmuji2(2)+a32*gmuji2(3)
+         geel_loc_ji=
+     &     +a22*gmuji2(1)
+     &     +a23*gmuji2(2)
+     &     +a32*gmuji2(3)
      &     +a33*gmuji2(4)
-         write(iout,*) "derivative over thataj-1"
-         write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3),
-     &   a33*gmuji2(4)
+c         write(iout,*) "derivative over thataj-1"
+c         write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3),
+c     &   a33*gmuji2(4)
          gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+
      &      geel_loc_ji*wel_loc
 #endif