working microcanonical for CA2+ protein

[unres4.git] / source / unres / geometry.f90

diff --git a/source/unres/geometry.f90 b/source/unres/geometry.f90
index 653e595..52c3b83 100644 (file)
--- a/source/unres/geometry.f90
+++ b/source/unres/geometry.f90
@@ -450,6 +450,7 @@
         alph(i)=alpha(nres+i,i,nres2+2)
         theta(i+1)=alpha(i-1,i,i+1)
         vbld(i)=dist(i-1,i)
+!        print *,i,vbld(i),"vbld(i)"
         vbld_inv(i)=1.0d0/vbld(i)
         vbld(nres+i)=dist(nres+i,i)
         if (itype(i,1).ne.10) then
@@ -2913,11 +2914,11 @@
        endif
       endif
       do i=1,nres-1
-       if (molnum(i).ne.1) cycle
+!       if (molnum(i).ne.1) cycle
 !in wham      do i=1,nres
         iti=itype(i,1)
-        if ((dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0).and.&
-       (iti.ne.ntyp1  .and. itype(i+1,1).ne.ntyp1)) then
+        if (((dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0).and.&
+       (iti.ne.ntyp1  .and. itype(i+1,1).ne.ntyp1)).and.molnum(i).eq.1) then
           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
 !test          stop
         endif
@@ -3246,7 +3247,8 @@
        do j=1,3
          gcart(j,1)=gcart(j,1)+gloc(1,icg)*dphi(j,1,4) &
            +gloc(nres-2,icg)*dtheta(j,1,3)      
-         if(itype(2,1).ne.10) then
+          if ((itype(2,1).ne.10).and.&
+          (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2)))) then
           gcart(j,1)=gcart(j,1)+gloc(ialph(2,1),icg)*dalpha(j,1,2)+ &
           gloc(ialph(2,1)+nside,icg)*domega(j,1,2)             
         endif
@@ -3273,11 +3275,15 @@
            gcart(j,3)=gcart(j,3)+gloc(1,icg)*dphi(j,3,4)+gloc(2,icg)* &
            dphi(j,2,5)+gloc(nres-1,icg)*dtheta(j,2,4)+gloc(nres,icg)* &
            dtheta(j,1,5)
-         if(itype(3,1).ne.10) then
+!         if(itype(3,1).ne.10) then
+          if ((itype(3,1).ne.10).and.&
+          (itype(3,molnum(3)).ne.ntyp1_molec(molnum(3)))) then
           gcart(j,3)=gcart(j,3)+gloc(ialph(3,1),icg)* &
           dalpha(j,2,3)+gloc(ialph(3,1)+nside,icg)*domega(j,2,3)
          endif
-        if(itype(4,1).ne.10) then
+!               if(itype(4,1).ne.10) then
+          if ((itype(4,1).ne.10).and.&
+          (itype(4,molnum(4)).ne.ntyp1_molec(molnum(4)))) then
           gcart(j,3)=gcart(j,3)+gloc(ialph(4,1),icg)* &
           dalpha(j,1,4)+gloc(ialph(4,1)+nside,icg)*domega(j,1,4)
          endif
@@ -3386,7 +3392,9 @@
        do j=1,3
          if((itype(2,1).ne.10).and. &
            (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2)))) then
-           if (itype(1,1).ne.10) gxcart(j,2)=gxcart(j,2)+ &
+           if ((itype(1,1).ne.10).and.&
+              ((itype(1,molnum(1)).ne.ntyp1_molec(molnum(1)))))&
+            gxcart(j,2)=gxcart(j,2)+ &
                                gloc_sc(3,0,icg)*dtauangle(j,3,3,3)
         if ((itype(3,1).ne.10).and.(nres.ge.3).and.(itype(3,molnum(3)).ne.ntyp1_molec(3))) &
          then
@@ -3395,6 +3403,8 @@
            gcart(j,2)=gcart(j,2)+gloc_sc(3,1,icg)*dtauangle(j,3,2,4)
           endif
           if (nres.gt.3) then
+!           if ((itype(1,1).ne.10).and.&
+!              ((itype(1,molnum(1)).ne.ntyp1_molec(molnum(1))))) &
            gxcart(j,2)=gxcart(j,2)-gloc_sc(1,1,icg)*dtauangle(j,1,1,4)
 !c                  the   - above is due to different vector direction
            gcart(j,2)=gcart(j,2)+gloc_sc(1,1,icg)*dtauangle(j,1,2,4)
@@ -3432,28 +3442,39 @@
           *dtauangle(j,1,2,i+2)
 !                   write(iout,*) "new",j,i,
 !     &  gcart(j,i),gloc_sc(1,i-1,icg),dtauangle(j,1,2,i+2)
-          if (itype(i-1,1).ne.10) then
+!          if (itype(i-1,1).ne.10) then
+          if ((itype(i-1,1).ne.10).and.&
+          (itype(i-1,molnum(i-1)).ne.ntyp1_molec(molnum(i-1)))) then
+
            gxcart(j,i)=gxcart(j,i)+gloc_sc(3,i-2,icg) &
       *dtauangle(j,3,3,i+1)
           endif
-          if (itype(i+1,1).ne.10) then
-           gxcart(j,i)=gxcart(j,i)-gloc_sc(3,i-1,icg) &
+!          if (itype(i+1,1).ne.10) then
+          if ((itype(i+1,1).ne.10).and.&
+          (itype(i+1,molnum(i+1)).ne.ntyp1_molec(molnum(i+1)))) then
+          gxcart(j,i)=gxcart(j,i)-gloc_sc(3,i-1,icg) &
       *dtauangle(j,3,1,i+2)
            gcart(j,i)=gcart(j,i)+gloc_sc(3,i-1,icg) &
       *dtauangle(j,3,2,i+2)
           endif
           endif
-          if (itype(i-1,1).ne.10) then
+!          if (itype(i-1,1).ne.10) then
+          if ((itype(i-1,1).ne.10).and.&
+          (itype(i-1,molnum(i-1)).ne.ntyp1_molec(molnum(i-1)))) then
            gcart(j,i)=gcart(j,i)+gloc_sc(1,i-2,icg)* &
            dtauangle(j,1,3,i+1)
           endif
-          if (itype(i+1,1).ne.10) then
+!          if (itype(i+1,1).ne.10) then
+          if ((itype(i+1,1).ne.10).and.&
+          (itype(i+1,molnum(i+1)).ne.ntyp1_molec(molnum(i+1)))) then
            gcart(j,i)=gcart(j,i)+gloc_sc(2,i-1,icg)* &
            dtauangle(j,2,2,i+2)
 !          write(iout,*) "numer",i,gloc_sc(2,i-1,icg),
 !     &    dtauangle(j,2,2,i+2)
           endif
-          if (itype(i+2,1).ne.10) then
+!          if (itype(i+2,1).ne.10) then
+          if ((itype(i+2,1).ne.10).and.&
+          (itype(i+2,molnum(i+2)).ne.ntyp1_molec(molnum(i+2)))) then
            gcart(j,i)=gcart(j,i)+gloc_sc(2,i,icg)* &
            dtauangle(j,2,1,i+3)
           endif
@@ -3469,8 +3490,10 @@
           *dtauangle(j,2,3,nres)
 !          write (iout,*) "gxcart(nres-1)", gloc_sc(2,nres-3,icg),
 !     &     dtauangle(j,2,3,nres), gxcart(j,nres-1)
-         if (itype(nres-2,1).ne.10) then
-        gxcart(j,nres-1)=gxcart(j,nres-1)+gloc_sc(3,nres-3,icg) &
+!         if (itype(nres-2,1).ne.10) then
+         if ((itype(nres-2,1).ne.10).and.&
+       (itype(nres-2,molnum(nres-2)).ne.ntyp1_molec(molnum(nres-2)))) then
+       gxcart(j,nres-1)=gxcart(j,nres-1)+gloc_sc(3,nres-3,icg) &
           *dtauangle(j,3,3,nres)
           endif
          if ((itype(nres,1).ne.10).and.&