X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2FREMD.f90;h=e2db47896619eb6a519ebd30fb2bccab65fb94e3;hb=705644e0cbb7678faefd6fe1bc436159d38ad85d;hp=1deb00096ba5c2446dc918e6026a38bd8379aca9;hpb=a01e4f4abfd7680ce30209afd6cab5d52532a55f;p=unres4.git

diff --git a/source/unres/REMD.f90 b/source/unres/REMD.f90
index 1deb000..e2db478 100644
--- a/source/unres/REMD.f90
+++ b/source/unres/REMD.f90
@@ -41,7 +41,7 @@
 !       include 'COMMON.CONTROL'
 !       include 'COMMON.MUCA'
 !      include 'COMMON.TIME1'
-       integer ::i,j,ind,itime,k
+       integer ::i,j,ind,itime,mnum
        real(kind=8) :: zapas(6*nres) !,muca_factor	!maxres6=6*maxres
        real(kind=8) :: rs(dimen),xsolv(dimen)
 #ifdef CHECK5DSOL
@@ -62,14 +62,16 @@
         enddo
         write (iout,*) "Potential forces sidechain"
         do i=nnt,nct
-          if (itype(i).ne.10 .and. itype(i).ne.ntyp1) &
+          mnum=molnum(i)
+          if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.mnum.eq.5)&
             write (iout,'(i5,3e15.5,5x,3e15.5)') i,(-gxcart(j,i),j=1,3)
         enddo
       endif
        do j=1,3
          ind=1
          do i=nnt,nct
-           if (itype(i).eq.10 .or. itype(i).eq.ntyp1)then
+           mnum=molnum(i)     
+           if (itype(i,1).eq.10 .or. itype(i,1).eq.ntyp1 .or. mnum.eq.5)then
 	    rs(ind)=-gcart(j,i)-gxcart(j,i)
             ind=ind+1
 	   else
@@ -122,7 +124,8 @@
 #endif
 	 ind=1
          do i=nnt,nct
-	  if (itype(i).eq.10 .or. itype(i).eq.ntyp1)then 
+          mnum=molnum(i)
+	  if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum))then 
 	   d_a(j,i)=xsolv(ind)
 	   ind=ind+1
 	  else
@@ -144,7 +147,9 @@
          d_a(j,0)=d_a(j,nnt)
        enddo
        do i=nnt,nct
-         if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then
+          mnum=molnum(i)
+!         if (itype(i,1).eq.10 .or. itype(i,1).eq.ntyp1) then
+          if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum))then
            do j=1,3
 	     d_a(j,i)=d_a(j,i+1)-d_a(j,i)
            enddo
@@ -188,7 +193,9 @@
       enddo
       if (lprn) write (iout,*) "Potential forces sidechain"
       do i=nnt,nct
-        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+          mnum=molnum(i)
+        if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
+          .and.mnum.ne.5) then
           if (lprn) write (iout,'(i5,3e15.5,5x,3e15.5)') &
              i,(-gxcart(j,i),j=1,3)
           do j=1,3
@@ -211,7 +218,10 @@
         enddo
       enddo
       do i=nnt,nct
-        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+       mnum=molnum(i)
+!        if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
+         if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
+          .and.mnum.ne.5) then
           do j=1,3
             ind=ind+1
             d_a(j,i+nres)=d_a_work(ind)
@@ -290,10 +300,10 @@
 !#endif
 !      include 'COMMON.IOUNITS'
 !      include 'COMMON.TIME1'
-      logical :: lprn = .true.
+      logical :: lprn = .false.
       logical :: osob
 #ifndef FIVEDIAG
-      real(kind=8),dimension(2*nres,2*nres) :: Bmat,matmult
+      real(kind=8),allocatable,dimension(:,:) :: Bmat,matmult
 #endif
       real(kind=8) :: dtdi
       real(kind=8),dimension(2*nres) :: massvec,sqreig	!(maxres2) maxres2=2*maxres
@@ -305,18 +315,22 @@
       real(kind=8),dimension(8*6*nres) :: work	!(8*maxres6)
       integer,dimension(6*nres) :: iwork	!(maxres6) maxres6=6*maxres
 !el      common /przechowalnia/ Gcopy,Ghalf
-      real(kind=8) :: coeff
+      real(kind=8) :: coeff,mscab
       integer :: i,j,ind,ind1,k,l,ii,jj,m,m1,ii1,iti,nres2,ierr,nind,mark
       real(kind=8) :: ip4
-      integer :: iz
+      integer :: iz,mnum
+      print *,"just entered"
       relfeh=1.0d-14
       nres2=2*nres
-
+      print *,"FIVEDIAG"
       write (iout,*) "before FIVEDIAG"
 #ifndef FIVEDIAG
       write (iout,*) "ALLOCATE"
+      print *,"ALLOCATE"
       if(.not.allocated(Gcopy)) allocate(Gcopy(nres2,nres2)) !(maxres2,maxres2)
       if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !mmaxres2=(maxres2*(maxres+1)/2)
+      if(.not.allocated(Bmat))  allocate(Bmat(nres2,nres2))
+      if(.not.allocated(matmult)) allocate(matmult(nres2,nres2))
 #endif
 !
 ! Set up the matrix of the (dC,dX)-->(C,X) transformation (A), the
@@ -329,6 +343,7 @@
       dimen3=dimen*3
       write (iout,*) "nnt",nnt," nct",nct," nside",nside
 #ifdef FIVEDIAG
+#ifdef CRYST_BOND
        ip4=ip/4
        DM(1)=mp/4+ip4
        if (iabs(itype(nnt)).eq.10) then
@@ -342,14 +357,14 @@
          nind=2
        endif	  
        do i=nnt+1,nct-1
-!         if (iabs(itype(i)).eq.ntyp1) cycle
+!         if (iabs(itype(i,1)).eq.ntyp1) cycle
          DM(ind)=2*ip4+mp/2
-         if (iabs(itype(i)).eq.10 .or. iabs(itype(i)).eq.ntyp1) then
-           if (iabs(itype(i)).eq.10) DM(ind)=DM(ind)+msc(10)
+         if (iabs(itype(i,1)).eq.10 .or. iabs(itype(i,1)).eq.ntyp1) then
+           if (iabs(itype(i,1)).eq.10) DM(ind)=DM(ind)+msc(10)
            ind=ind+1
          else
-           DM(ind)=DM(ind)+isc(iabs(itype(i)))
-           DM(ind+1)=msc(iabs(itype(i)))+isc(iabs(itype(i)))
+           DM(ind)=DM(ind)+isc(iabs(itype(i,1)))
+           DM(ind+1)=msc(iabs(itype(i,1)))+isc(iabs(itype(i,1)))
            ind=ind+2
          endif 
        enddo
@@ -368,8 +383,10 @@
 	
         ind=1
         do i=nnt,nct
-	if (iabs(itype(i)).ne.10 .and.iabs(itype((i))).ne.ntyp1) then
-          DU1(ind)=-isc(iabs(itype(i)))
+        mnum=molnum(i)
+	if (iabs(itype(i,1)).ne.10 .and.iabs(itype((i))).ne.ntyp1 &
+          .and.mnum.eq.5) then
+          DU1(ind)=-isc(iabs(itype(i,1)))
           DU1(ind+1)=0.0d0
 	  ind=ind+2
 	else
@@ -380,9 +397,11 @@
 
        ind=1
        do i=nnt,nct-1
-!        if (iabs(itype(i)).eq.ntyp1) cycle
-        write (iout,*) "i",i," itype",itype(i),ntyp1
-	if (iabs(itype(i)).ne.10 .and. iabs(itype(i)).ne.ntyp1) then
+       mnum=molnum(i)
+!        if (iabs(itype(i,1)).eq.ntyp1) cycle
+        write (iout,*) "i",i," itype",itype(i,1),ntyp1
+	if (iabs(itype(i,1)).ne.10 .and. &
+          iabs(itype(i,mnum)).ne.ntyp1_molec(mnum) .and. mnum.ne.5) then
          DU2(ind)=mp/4-ip4
          DU2(ind+1)=0.0d0
          ind=ind+2
@@ -392,6 +411,75 @@
 	 ind=ind+1
 	endif
        enddo
+#else
+       ip4=ip(1)/4
+       DM(1)=mp(1)/4+ip4
+       if (iabs(itype(nnt,1)).eq.10) then
+         DM(1)=DM(1)+msc(10,1)
+         ind=2
+         nind=1
+       else
+         DM(1)=DM(1)+isc(iabs(itype(nnt,1)),1)
+         DM(2)=msc(iabs(itype(nnt,1)),1)+isc(iabs(itype(nnt,1)),1)
+         ind=3
+         nind=2
+       endif     
+       do i=nnt+1,nct-1
+!         if (iabs(itype(i,1)).eq.ntyp1) cycle
+         DM(ind)=2*ip4+mp(1)/2
+         if (iabs(itype(i,1)).eq.10 .or. iabs(itype(i,1)).eq.ntyp1) then
+           if (iabs(itype(i,1)).eq.10) DM(ind)=DM(ind)+msc(10,1)
+           ind=ind+1
+         else
+           DM(ind)=DM(ind)+isc(iabs(itype(i,1)),1)
+           DM(ind+1)=msc(iabs(itype(i,1)),1)+isc(iabs(itype(i,1)),1)
+           ind=ind+2
+         endif 
+       enddo
+       if (nct.gt.nnt) then
+       DM(ind)=ip4+mp(1)/4
+       if (iabs(itype(nct,1)).eq.10) then
+         DM(ind)=DM(ind)+msc(10,1)
+         nind=ind
+       else
+         DM(ind)=DM(ind)+isc(iabs(itype(nct,1)),1)
+         DM(ind+1)=msc(iabs(itype(nct,1)),1)+isc(iabs(itype(nct,1)),1)
+         nind=ind+1
+       endif
+        endif
+       
+       
+        ind=1
+        do i=nnt,nct
+        mnum=molnum(i)
+       if (iabs(itype(i,1)).ne.10 .and.iabs(itype(i,mnum)).ne.ntyp1 &
+          .and.mnum.eq.5) then
+          DU1(ind)=-isc(iabs(itype(i,1)),1)
+          DU1(ind+1)=0.0d0
+         ind=ind+2
+       else
+         DU1(ind)=mp(1)/4-ip4
+         ind=ind+1
+       endif
+       enddo
+
+       ind=1
+       do i=nnt,nct-1
+       mnum=molnum(i)
+!        if (iabs(itype(i,1)).eq.ntyp1) cycle
+        write (iout,*) "i",i," itype",itype(i,1),ntyp1
+       if (iabs(itype(i,1)).ne.10 .and. &
+          iabs(itype(i,mnum)).ne.ntyp1_molec(mnum) .and. mnum.ne.5) then
+         DU2(ind)=mp(1)/4-ip4
+         DU2(ind+1)=0.0d0
+         ind=ind+2
+       else
+         DU2(ind)=0.0d0
+        DU2(ind+1)=0.0d0
+        ind=ind+1
+       endif
+       enddo
+#endif
        DMorig=DM
        DU1orig=DU1
        DU2orig=DU2
@@ -499,12 +587,17 @@
 !  Diagonal elements of the dC part of A and the respective friction coefficients
       ind=1
       ind1=0
+      print *,"TUTUTUT?!",nnt,nct-1
       do i=nnt,nct-1
+        mnum=molnum(i)
         ind=ind+1
         ind1=ind1+1
-        coeff=0.25d0*IP
-        massvec(ind1)=mp
+        coeff=0.25d0*IP(mnum)
+        if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
+        print *,i,coeff,mp(mnum)
+        massvec(ind1)=mp(mnum)
         Gmat(ind,ind)=coeff
+        print *,"i",mp(mnum)
         A(ind1,ind)=0.5d0
       enddo
       
@@ -520,24 +613,33 @@
       m1=nct-nnt+1
       ind=0
       ind1=0
-      msc(ntyp1)=1.0d0
+      do i=1,2
+      msc(ntyp1_molec(i),i)=1.0d0
+      enddo
       do i=nnt,nct
         ind=ind+1
         ii = ind+m
-        iti=itype(i)
-        massvec(ii)=msc(iabs(iti))
-        if (iti.ne.10 .and. iti.ne.ntyp1) then
+        mnum=molnum(i)
+        iti=itype(i,mnum)
+        if (mnum.eq.5) then
+        mscab=0.0
+        else
+        mscab=msc(iabs(iti),mnum)
+        endif
+        massvec(ii)=mscab
+
+        if (iti.ne.10 .and. iti.ne.ntyp1_molec(mnum) .and. mnum.ne.5) then
           ind1=ind1+1
           ii1= ind1+m1
           A(ii,ii1)=1.0d0
-          Gmat(ii1,ii1)=ISC(iabs(iti))
+          Gmat(ii1,ii1)=ISC(iabs(iti),mnum)
         endif
       enddo
 !  Off-diagonal elements of the dX part of A
       ind=0
       k=nct-nnt
       do i=nnt,nct
-        iti=itype(i)
+        iti=itype(i,1)
         ind=ind+1
         do j=nnt,i
           ii = ind
@@ -588,7 +690,9 @@
       Bmat(1,1)=1.0d0
       j=2
       do i=nnt,nct
-        if (itype(i).eq.10) then
+        mnum=molnum(i)
+        if ((itype(i,1).eq.10).or.(itype(i,mnum).eq.ntyp1_molec(mnum))&
+          .or.(mnum.eq.5)) then
           Bmat(i-nnt+2,j-1)=-1
           Bmat(i-nnt+2,j)=1
           j=j+1
@@ -611,15 +715,20 @@
         call matout(dimen,dimen,nres2,nres2,Bmat)
       endif
       Gcopy=0.0d0
+      write(iout,*) "before Gcopy",dimen,nres*2
+#ifdef TESTFIVEDIAG
       do i=1,dimen
         do j=1,dimen
+!            write (iout,*) "myij",i,j
             do k=1,dimen
                do l=1,dimen
+!                 write(iout,*) "i,j,k,l",i,j,k,l
                  Gcopy(i,j)=Gcopy(i,j)+Bmat(k,i)*Gmat(k,l)*Bmat(l,j)
                enddo
             enddo
         enddo
       enddo
+#endif
       if (lprn) then
         write (iout,'(//a)') "Gmat-transformed"
         call matout(dimen,dimen,nres2,nres2,Gcopy)