working microcanonical
[unres4.git] / source / unres / geometry.f90
index e0b37d3..52c3b83 100644 (file)
         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
        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
 !   calculating dE/ddc1  
 !el local variables
        integer :: j,i
+!       print *,"gloc",gloc(:,:)
+!       print *, "gcart",gcart(:,:)
        if (nres.lt.3) go to 18
        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
            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
         enddo
 !   The side-chain vector derivatives
         do i=2,nres-1
-         if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then   
+         if(itype(i,1).ne.10 .and.  &
+           itype(i,molnum(i)).ne.ntyp1_molec(molnum(i))) then  
             do j=1,3   
               gxcart(j,i)=gxcart(j,i)+gloc(ialph(i,1),icg)*dalpha(j,3,i) &
               +gloc(ialph(i,1)+nside,icg)*domega(j,3,i)
 !       enddo
        if (nres.lt.2) return
        if ((nres.lt.3).and.(itype(1,1).eq.10)) return
-       if ((itype(1,1).ne.10).and.(itype(1,1).ne.ntyp1)) then
+       if ((itype(1,1).ne.10).and. &
+        (itype(1,molnum(1)).ne.ntyp1_molec(molnum(1)))) then
         do j=1,3
 !c Derviative was calculated for oposite vector of side chain therefore
 ! there is "-" sign before gloc_sc
            dtauangle(j,1,1,3)
          gcart(j,1)=gcart(j,1)+gloc_sc(1,0,icg)* &
            dtauangle(j,1,2,3)
-          if ((itype(2,1).ne.10).and.(itype(2,1).ne.ntyp1)) then
+          if ((itype(2,1).ne.10).and. &
+        (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2)))) then
          gxcart(j,1)= gxcart(j,1) &
                      -gloc_sc(3,0,icg)*dtauangle(j,3,1,3)
          gcart(j,1)=gcart(j,1)+gloc_sc(3,0,icg)* &
           endif
        enddo
        endif
-         if ((nres.ge.3).and.(itype(3,1).ne.10).and.(itype(3,1).ne.ntyp1)) &
+         if ((nres.ge.3).and.(itype(3,molnum(3)).ne.10).and.&
+         (itype(3,molnum(3)).ne.ntyp1_molec(molnum(3)))) &
       then
          do j=1,3
          gcart(j,1)=gcart(j,1)+gloc_sc(2,1,icg)*dtauangle(j,2,1,4)
 
 !     Calculating the remainder of dE/ddc2
        do j=1,3
-         if((itype(2,1).ne.10).and.(itype(2,1).ne.ntyp1)) then
-           if (itype(1,1).ne.10) gxcart(j,2)=gxcart(j,2)+ &
+         if((itype(2,1).ne.10).and. &
+           (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2)))) 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(3,0,icg)*dtauangle(j,3,3,3)
-        if ((itype(3,1).ne.10).and.(nres.ge.3).and.(itype(3,1).ne.ntyp1)) &
+        if ((itype(3,1).ne.10).and.(nres.ge.3).and.(itype(3,molnum(3)).ne.ntyp1_molec(3))) &
          then
            gxcart(j,2)=gxcart(j,2)-gloc_sc(3,1,icg)*dtauangle(j,3,1,4)
 !c                  the   - above is due to different vector direction
            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)
 !           write(iout,*) gloc_sc(1,1,icg),dtauangle(j,1,1,4),"gx"
           endif
          endif
-         if ((itype(1,1).ne.10).and.(itype(1,1).ne.ntyp1)) then
+         if ((itype(1,1).ne.10).and.&
+         (itype(1,molnum(1)).ne.ntyp1_molec(molnum(1)))) then
           gcart(j,2)=gcart(j,2)+gloc_sc(1,0,icg)*dtauangle(j,1,3,3)
 !           write(iout,*)  gloc_sc(1,0,icg),dtauangle(j,1,3,3)
         endif
         do i=3,nres-2
          do j=1,3
 !          write(iout,*) "before", gcart(j,i)
-          if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then
+          if ((itype(i,1).ne.10).and.&
+          (itype(i,molnum(i)).ne.ntyp1_molec(molnum(i)))) then
           gxcart(j,i)=gxcart(j,i)+gloc_sc(2,i-2,icg) &
           *dtauangle(j,2,3,i+1) &
           -gloc_sc(1,i-1,icg)*dtauangle(j,1,1,i+2)
           *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
 !  Setting dE/ddnres-1       
       if(nres.ge.4) then
          do j=1,3
-         if ((itype(nres-1,1).ne.10).and.(itype(nres-1,1).ne.ntyp1)) then
+         if ((itype(nres-1,1).ne.10).and.&
+       (itype(nres-1,molnum(nres-1)).ne.ntyp1_molec(molnum(nres-1)))) then
          gxcart(j,nres-1)=gxcart(j,nres-1)+gloc_sc(2,nres-3,icg) &
           *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.(itype(nres,1).ne.ntyp1)) then
+         if ((itype(nres,1).ne.10).and.&
+         (itype(nres,molnum(nres)).ne.ntyp1_molec(molnum(nres)))) then
         gxcart(j,nres-1)=gxcart(j,nres-1)-gloc_sc(3,nres-2,icg) &
           *dtauangle(j,3,1,nres+1)
         gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(3,nres-2,icg) &
           *dtauangle(j,3,2,nres+1)
           endif
          endif
-         if ((itype(nres-2,1).ne.10).and.(itype(nres-2,1).ne.ntyp1)) then
+         if ((itype(nres-2,1).ne.10).and.&
+         (itype(nres-2,molnum(nres-2)).ne.ntyp1_molec(molnum(nres-2)))) then
             gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(1,nres-3,icg)* &
          dtauangle(j,1,3,nres)
          endif
-          if ((itype(nres,1).ne.10).and.(itype(nres,1).ne.ntyp1)) then
+          if ((itype(nres,1).ne.10).and.(itype(nres,molnum(nres)).ne.ntyp1_molec(molnum(nres)))) then
             gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(2,nres-2,icg)* &
            dtauangle(j,2,2,nres+1)
 !           write (iout,*) "gcart(nres-1)", gloc_sc(2,nres-2,icg),
       endif
 !  Settind dE/ddnres       
        if ((nres.ge.3).and.(itype(nres,1).ne.10).and. &
-          (itype(nres,1).ne.ntyp1))then
+          (itype(nres,molnum(nres)).ne.ntyp1_molec(molnum(nres))))then
        do j=1,3
         gxcart(j,nres)=gxcart(j,nres)+gloc_sc(3,nres-2,icg) &
        *dtauangle(j,3,3,nres+1)+gloc_sc(2,nres-2,icg) &
         enddo
        endif
 !   The side-chain vector derivatives
+!       print *,"gcart",gcart(:,:)
       return
       end subroutine int_to_cart
 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)