changes in ions
[unres4.git] / source / unres / geometry.F90
index 9641e8c..ae86ead 100644 (file)
       chiom1=chi1*om1
       chiom2=chi2*om2
       facsig=om1*chiom1+om2*chiom2-2.0D0*om1om2*chiom12
+!      print *,"TUT?",om1*chiom1,facsig,om1,om2,om12
       sigsq=1.0D0-facsig*faceps1_inv
       sigsq_om1=(chiom1-chiom12*om2)*faceps1_inv
       sigsq_om2=(chiom2-chiom12*om1)*faceps1_inv
        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)      
+!         write(iout,*) "pierwszy gcart", gcart(j,2)
           if ((itype(2,1).ne.10).and.&
-          (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2)))) then
+          (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2)).and.(molnum(2).ne.5))) 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
        do j=1,3
          gcart(j,2)=gcart(j,2)+gloc(1,icg)*dphi(j,2,4)+ &
         gloc(nres-2,icg)*dtheta(j,2,3)+gloc(nres-1,icg)*dtheta(j,1,4)
-        if(itype(2,1).ne.10) then
+!         write(iout,*) "drugi gcart", gcart(j,2)
+        if((itype(2,1).ne.10).and.(molnum(2).ne.5)) then
           gcart(j,2)=gcart(j,2)+gloc(ialph(2,1),icg)*dalpha(j,2,2)+ &
           gloc(ialph(2,1)+nside,icg)*domega(j,2,2)
         endif
 !   The side-chain vector derivatives
         do i=2,nres-1
          if(itype(i,1).ne.10 .and.  &
-           itype(i,molnum(i)).ne.ntyp1_molec(molnum(i))) then  
+           itype(i,molnum(i)).ne.ntyp1_molec(molnum(i)).and.&
+             (molnum(i).ne.5)) 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)
 ! INTERTYP=3 SC...Ca...Ca...SC
 !   calculating dE/ddc1      
   18   continue
+!         write(iout,*) "przed sccor gcart", gcart(1,2)
+
 !       do i=1,nres
 !       gloc(i,icg)=0.0D0
 !          write (iout,*) "poczotkoawy",i,gloc_sc(1,i,icg)
        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,molnum(1)).ne.ntyp1_molec(molnum(1)))) then
+        (itype(1,molnum(1)).ne.ntyp1_molec(molnum(1))).and.(molnum(1).ne.5)) then
         do j=1,3
 !c Derviative was calculated for oposite vector of side chain therefore
 ! there is "-" sign before gloc_sc
          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,molnum(2)).ne.ntyp1_molec(molnum(2)))) then
+        (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2))).and.(molnum(2).ne.5)) 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)* &
 !   ommited 
 !     &     +gloc_sc(intertyp,nres-2,icg)*dtheta(j,1,3)  
 
+!         write(iout,*) "przed dE/ddc2 gcart", gcart(1,2)
+
 !     Calculating the remainder of dE/ddc2
        do j=1,3
          if((itype(2,1).ne.10).and. &
-           (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2)))) then
+           (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2)).and.(molnum(2).ne.5))) then
            if ((itype(1,1).ne.10).and.&
-              ((itype(1,molnum(1)).ne.ntyp1_molec(molnum(1)))))&
+              ((itype(1,molnum(1)).ne.ntyp1_molec(molnum(1)))).and.(molnum(1).ne.5))&
             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))) &
+        if ((itype(3,1).ne.10).and.(nres.ge.3).and.(itype(3,molnum(3)).ne.ntyp1_molec(3)).and.molnum(3).ne.5) &
          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
            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,2,4),"gcart"
+!          write(iout,*) gloc_sc(1,1,icg),dtauangle(j,1,2,4),"gcart",gcart(j,2)
 !           write(iout,*) gloc_sc(1,1,icg),dtauangle(j,1,1,4),"gx"
           endif
          endif
          if ((itype(1,1).ne.10).and.&
-         (itype(1,molnum(1)).ne.ntyp1_molec(molnum(1)))) then
+         (itype(1,molnum(1)).ne.ntyp1_molec(molnum(1))).and.(molnum(1).ne.5)) 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
-         if ((itype(3,1).ne.10).and.(nres.ge.3)) then
+         if ((itype(3,1).ne.10).and.(nres.ge.3).and.(molnum(3).ne.5)) then
           gcart(j,2)=gcart(j,2)+gloc_sc(2,1,icg)*dtauangle(j,2,2,4)
 !           write(iout,*) gloc_sc(2,1,icg),dtauangle(j,2,2,4)
          endif
-         if ((itype(4,1).ne.10).and.(nres.ge.4)) then
+         if ((itype(4,1).ne.10).and.(nres.ge.4).and.(molnum(4).ne.5)) then
           gcart(j,2)=gcart(j,2)+gloc_sc(2,2,icg)*dtauangle(j,2,1,5)
 !           write(iout,*) gloc_sc(2,2,icg),dtauangle(j,2,1,5)
          endif
 
 !      write(iout,*) gcart(j,2),itype(2,1),itype(1,1),itype(3,1), "gcart2"
        enddo
+!         write(iout,*) "po dE/ddc2 gcart", gcart(1,2)
+
 !    If there are more than five residues
       if(nres.ge.5) then                        
         do i=3,nres-2
          do j=1,3
 !          write(iout,*) "before", gcart(j,i)
           if ((itype(i,1).ne.10).and.&
-          (itype(i,molnum(i)).ne.ntyp1_molec(molnum(i)))) then
+          (itype(i,molnum(i)).ne.ntyp1_molec(molnum(i))).and.(molnum(i).ne.5)) 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)
 !     &  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).and.&
-          (itype(i-1,molnum(i-1)).ne.ntyp1_molec(molnum(i-1)))) then
+          (itype(i-1,molnum(i-1)).ne.ntyp1_molec(molnum(i-1))).and.(molnum(i-1).eq.5)) then
 
            gxcart(j,i)=gxcart(j,i)+gloc_sc(3,i-2,icg) &
       *dtauangle(j,3,3,i+1)
           endif
 !          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
+          (itype(i+2,molnum(i+2)).ne.ntyp1_molec(molnum(i+2))).and.(molnum(i+2).ne.5)) then
            gcart(j,i)=gcart(j,i)+gloc_sc(2,i,icg)* &
            dtauangle(j,2,1,i+3)
           endif
       if(nres.ge.4) then
          do j=1,3
          if ((itype(nres-1,1).ne.10).and.&
-       (itype(nres-1,molnum(nres-1)).ne.ntyp1_molec(molnum(nres-1)))) then
+       (itype(nres-1,molnum(nres-1)).ne.ntyp1_molec(molnum(nres-1))).and.(molnum(nres-1).ne.5)) 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
          if ((itype(nres-2,1).ne.10).and.&
-       (itype(nres-2,molnum(nres-2)).ne.ntyp1_molec(molnum(nres-2)))) then
+       (itype(nres-2,molnum(nres-2)).ne.ntyp1_molec(molnum(nres-2))).and.(molnum(nres-2).ne.5)) 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,molnum(nres)).ne.ntyp1_molec(molnum(nres)))) then
+         (itype(nres,molnum(nres)).ne.ntyp1_molec(molnum(nres))).and.(molnum(nres).ne.5)) 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) &
           endif
          endif
          if ((itype(nres-2,1).ne.10).and.&
-         (itype(nres-2,molnum(nres-2)).ne.ntyp1_molec(molnum(nres-2)))) then
+         (itype(nres-2,molnum(nres-2)).ne.ntyp1_molec(molnum(nres-2))).and.(molnum(nres-2).ne.5)) 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,molnum(nres)).ne.ntyp1_molec(molnum(nres)))) then
+          if ((itype(nres,1).ne.10).and.(itype(nres,molnum(nres)).ne.ntyp1_molec(molnum(nres))).and.(molnum(nres).ne.5)) 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,molnum(nres)).ne.ntyp1_molec(molnum(nres))))then
+          (itype(nres,molnum(nres)).ne.ntyp1_molec(molnum(nres))).and.(molnum(nres).ne.5))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) &
        *dtauangle(j,2,3,nres+1)
         enddo
        endif
+!       write(iout,*) "final gcart",gcart(1,2)
 !   The side-chain vector derivatives
 !       print *,"gcart",gcart(:,:)
       return
         do i=1,nres
          if (molnum(i).eq.5) then
           c(1,i)=dmod(c(1,i),boxxsize)
+          if (c(1,i).lt.0) c(1,i)=c(1,i)+boxxsize
           c(2,i)=dmod(c(2,i),boxysize)
+          if (c(2,i).lt.0) c(2,i)=c(2,i)+boxysize
           c(3,i)=dmod(c(3,i),boxzsize)
+          if (c(3,i).lt.0) c(3,i)=c(3,i)+boxzsize
           c(1,i+nres)=dmod(c(1,i+nres),boxxsize)
           c(2,i+nres)=dmod(c(2,i+nres),boxysize)
           c(3,i+nres)=dmod(c(3,i+nres),boxzsize)