X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2FREMD.f90;h=1e294adfc1047114a57826de290044372ce756cf;hb=8ca97b16fe25b7053f258263899ba030572cc58f;hp=fe970be5ab13e7377f8fe65f8152406f2ff62082;hpb=d0aebe86703877377e5e16d63927e242b0e23469;p=unres4.git

diff --git a/source/unres/REMD.f90 b/source/unres/REMD.f90
index fe970be..1e294ad 100644
--- a/source/unres/REMD.f90
+++ b/source/unres/REMD.f90
@@ -41,13 +41,13 @@
 !       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
        real(kind=8) :: rscheck(dimen),rsold(dimen)
 #endif
-       logical :: lprn = .false.
+       logical :: lprn = .true.
 !el       common /cipiszcze/ itime
        itime = itt_comm
 
@@ -62,14 +62,14 @@
         enddo
         write (iout,*) "Potential forces sidechain"
         do i=nnt,nct
-          if (itype(i).ne.10 .and. itype(i).ne.ntyp1) &
+          if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) &
             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
+           if (itype(i,1).eq.10 .or. itype(i,1).eq.ntyp1)then
 	    rs(ind)=-gcart(j,i)-gxcart(j,i)
             ind=ind+1
 	   else
@@ -122,7 +122,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 +145,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 +191,8 @@
       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)) then
           if (lprn) write (iout,'(i5,3e15.5,5x,3e15.5)') &
              i,(-gxcart(j,i),j=1,3)
           do j=1,3
@@ -211,7 +215,9 @@
         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)) then
           do j=1,3
             ind=ind+1
             d_a(j,i+nres)=d_a_work(ind)
@@ -290,7 +296,7 @@
 !#endif
 !      include 'COMMON.IOUNITS'
 !      include 'COMMON.TIME1'
-      logical :: lprn = .false.
+      logical :: lprn = .true.
       logical :: osob
 #ifndef FIVEDIAG
       real(kind=8),dimension(2*nres,2*nres) :: Bmat,matmult
@@ -308,7 +314,7 @@
       real(kind=8) :: coeff
       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
       relfeh=1.0d-14
       nres2=2*nres
 
@@ -342,14 +348,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 +374,8 @@
 	
         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)))
+	if (iabs(itype(i,1)).ne.10 .and.iabs(itype((i))).ne.ntyp1) then
+          DU1(ind)=-isc(iabs(itype(i,1)))
           DU1(ind+1)=0.0d0
 	  ind=ind+2
 	else
@@ -380,9 +386,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
@@ -499,12 +507,15 @@
 !  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)
+        massvec(ind1)=mp(mnum)
         Gmat(ind,ind)=coeff
+        print *,"i",mp(mnum)
         A(ind1,ind)=0.5d0
       enddo
       
@@ -520,24 +531,27 @@
       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)
+        massvec(ii)=msc(iabs(iti),mnum)
+        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 +602,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