Merge branch 'UCGM' of mmka.chem.univ.gda.pl:unres4 into UCGM

[unres4.git] / source / unres / REMD.f90

diff --git a/source/unres/REMD.f90 b/source/unres/REMD.f90
index 89fddc9..3d11d5c 100644 (file)
--- a/source/unres/REMD.f90
+++ b/source/unres/REMD.f90
@@ -343,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
@@ -410,6 +411,77 @@
         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
+       if (iabs(itype(i,1)).ne.10 .and. &
+          iabs(itype(i,mnum)).ne.ntyp1_molec(mnum) .and. mnum.ne.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
@@ -517,7 +589,7 @@
 !  Diagonal elements of the dC part of A and the respective friction coefficients
       ind=1
       ind1=0
-      print *,"TUTUTUT?!",nnt,nct-1
+!      print *,"TUTUTUT?!",nnt,nct-1
       do i=nnt,nct-1
         mnum=molnum(i)
         ind=ind+1
@@ -527,7 +599,7 @@
         print *,i,coeff,mp(mnum)
         massvec(ind1)=mp(mnum)
         Gmat(ind,ind)=coeff
-        print *,"i",mp(mnum)
+!        print *,"i",mp(mnum)
         A(ind1,ind)=0.5d0
       enddo
       
@@ -738,6 +810,7 @@
       endif
       deallocate(Gcopy)
 #endif
+!write(iout,*) "end setup_MD_matr"
       return
       end subroutine setup_MD_matrices
 !-----------------------------------------------------------------------------