poprawki w ANG w src_MD-M

diff --git a/source/unres/src_MD-M/arcos.f b/source/unres/src_MD-M/arcos.f
index f054118..1e355ec 100644 (file)
--- a/source/unres/src_MD-M/arcos.f
+++ b/source/unres/src_MD-M/arcos.f
@@ -1,4 +1,4 @@
-      FUNCTION ARCOS(X)
+      double precision FUNCTION ARCOS(X)
       implicit real*8 (a-h,o-z)
       include 'COMMON.GEO'
       IF (DABS(X).LT.1.0D0) GOTO 1

diff --git a/source/unres/src_MD-M/checkder_p.F b/source/unres/src_MD-M/checkder_p.F
index 0539e48..b017e1c 100644 (file)
--- a/source/unres/src_MD-M/checkder_p.F
+++ b/source/unres/src_MD-M/checkder_p.F
@@ -281,7 +281,7 @@ C Check the gradient of the energy in Cartesian coordinates.
 c      call intcartderiv
 c      call checkintcartgrad
       call zerograd
-      aincr=1.0D-5
+      aincr=10.0D-6
       write(iout,*) 'Calling CHECK_ECARTINT.'
       nf=0
       icall=0
@@ -362,6 +362,7 @@ c      call checkintcartgrad
       endif
       write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
       do i=0,nres
+        print *,i
         do j=1,3
          xx(j)=c(j,i+nres)
          ddc(j)=dc(j,i) 

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 fa6c1de..e7199c4 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -4574,7 +4574,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 +4588,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 +4605,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 +4626,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 +4676,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 +4697,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 +4742,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 +4751,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