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

diff --git a/source/unres/REMD.f90 b/source/unres/REMD.f90
index 8a049a7..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
 
@@ -122,7 +122,8 @@
 #endif
 	 ind=1
          do i=nnt,nct
-	  if (itype(i,1).eq.10 .or. itype(i,1).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,1).eq.10 .or. itype(i,1).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,1).ne.10 .and. itype(i,1).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,1).ne.10 .and. itype(i,1).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
 
@@ -380,9 +386,11 @@
 
        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,1)).ne.ntyp1) then
+	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(1)
-        massvec(ind1)=mp(1)
+        coeff=0.25d0*IP(mnum)
+        massvec(ind1)=mp(mnum)
         Gmat(ind,ind)=coeff
+        print *,"i",mp(mnum)
         A(ind1,ind)=0.5d0
       enddo
       
@@ -520,17 +531,20 @@
       m1=nct-nnt+1
       ind=0
       ind1=0
-      msc(ntyp1,1)=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,1)
-        massvec(ii)=msc(iabs(iti),1)
-        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),1)
+          Gmat(ii1,ii1)=ISC(iabs(iti),mnum)
         endif
       enddo
 !  Off-diagonal elements of the dX part of A
@@ -588,7 +602,9 @@
       Bmat(1,1)=1.0d0
       j=2
       do i=nnt,nct
-        if (itype(i,1).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