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

diff --git a/source/unres/REMD.f90 b/source/unres/REMD.f90
index e1cf84c..1e294ad 100644
--- a/source/unres/REMD.f90
+++ b/source/unres/REMD.f90
@@ -25,7 +25,7 @@
       use comm_cipiszcze
       use energy_data
       use geometry_data, only: nres
-      use control_data 
+      use control_data    !el, only: mucadyn,lmuca
 #ifdef MPI
        include 'mpif.h'
       real(kind=8) :: time00
@@ -41,16 +41,139 @@
 !       include 'COMMON.CONTROL'
 !       include 'COMMON.MUCA'
 !      include 'COMMON.TIME1'
-       
-       integer :: i,j,ind,itime
+       integer ::i,j,ind,itime,mnum
        real(kind=8) :: zapas(6*nres) !,muca_factor	!maxres6=6*maxres
-       logical :: lprn = .false.
+       real(kind=8) :: rs(dimen),xsolv(dimen)
+#ifdef CHECK5DSOL
+       real(kind=8) :: rscheck(dimen),rsold(dimen)
+#endif
+       logical :: lprn = .true.
 !el       common /cipiszcze/ itime
        itime = itt_comm
 
 #ifdef TIMING
        time00=MPI_Wtime()
 #endif
+#ifdef FIVEDIAG
+      if (lprn) then
+        write (iout,*) "Potential forces backbone"
+        do i=nnt,nct
+          write (iout,'(i5,3e15.5,5x,3e15.5)')i,(-gcart(j,i),j=1,3)
+        enddo
+        write (iout,*) "Potential forces sidechain"
+        do i=nnt,nct
+          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,1).eq.10 .or. itype(i,1).eq.ntyp1)then
+	    rs(ind)=-gcart(j,i)-gxcart(j,i)
+            ind=ind+1
+	   else
+	    rs(ind)=-gcart(j,i)
+	    rs(ind+1)=-gxcart(j,i)
+	    ind=ind+2
+	   end if 
+         enddo
+#ifdef CHECK5DSOL
+         rsold=rs 
+#endif
+         if (lprn) then
+           write(iout,*) "RHS of the 5-diag equations system"
+           do i=1,dimen
+             write(iout,*) i,rs(i)
+           enddo
+         endif
+	 
+         call FDISYS (dimen,DM,DU1,DU2,rs,xsolv)
+         if (lprn) then
+	   write (iout,*) "Solution of the 5-diagonal equations system"
+           do i=1,dimen
+             write (iout,'(i5,f10.5)') i,xsolv(i)
+	   enddo
+         endif
+#ifdef CHECK5DSOL
+! Check the solution
+         rscheck(1)=DMorig(1)*xsolv(1)+DU1orig(1)*xsolv(2)+&
+           DU2orig(1)*xsolv(3) 
+         rscheck(2)=DU1orig(1)*xsolv(1)+DMorig(2)*xsolv(2)+&
+           DU1orig(2)*xsolv(3)+DU2orig(2)*xsolv(4)
+         
+         do i=3,dimen-2
+          rscheck(i)=DU2orig(i-2)*xsolv(i-2)+DU1orig(i-1)*&
+            xsolv(i-1)+DMorig(i)*xsolv(i)+DU1orig(i)*&
+             xsolv(i+1)+DU2orig(i)*xsolv(i+2)
+         enddo
+       rscheck(dimen-1)=DU2orig(dimen-3)*xsolv(dimen-3)+&
+         DU1orig(dimen-2)*xsolv(dimen-2)+DMorig(dimen-1)*&
+          xsolv(dimen-1)+DU1orig(dimen-1)*&
+           xsolv(dimen)
+       rscheck(dimen)=DU2orig(dimen-2)*xsolv(dimen-2)+DU1orig(dimen-1)*&
+          xsolv(dimen-1)+DMorig(dimen)*xsolv(dimen)
+         
+         do i=1,dimen
+            write(iout,*) "i",i,"rsold",rsold(i),"rscheck",rscheck(i),&
+            "ratio",rscheck(i)/rsold(i)
+           enddo
+! end check
+#endif
+	 ind=1
+         do i=nnt,nct
+          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
+	   d_a(j,i)=xsolv(ind)
+	   d_a(j,i+nres)=xsolv(ind+1)
+	  ind=ind+2
+	  end if 
+	 enddo
+       enddo
+       if (lprn) then
+       do i=nnt,nct
+        write (iout,'(i3,3f10.5)') i,(d_a(j,i),j=1,3)
+       enddo
+       do i=nnt,nct
+        write (iout,'(i3,3f10.5)') i,(d_a(j,i+nres),j=1,3)
+       enddo
+       endif
+       do j=1,3
+         d_a(j,0)=d_a(j,nnt)
+       enddo
+       do i=nnt,nct
+          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
+	 else 
+           do j=1,3
+	     d_a(j,i+nres)=d_a(j,i+nres)-d_a(j,i)
+	     d_a(j,i)=d_a(j,i+1)-d_a(j,i)
+           enddo
+       
+
+	 end if
+       enddo
+
+      if (lprn) then
+        write(iout,*) 'acceleration 3D FIVEDIAG'
+        write (iout,'(i3,3f10.5,3x,3f10.5)') 0,(d_a(j,0),j=1,3)
+        do i=nnt,nct-1
+         write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3)
+        enddo
+        do i=nnt,nct
+         write (iout,'(i3,3f10.5,3x,3f10.5)') &
+           i+nres,(d_a(j,i+nres),j=1,3)
+        enddo
+      endif
+#else
+! Old procedure
        do j=1,3
          zapas(j)=-gcart(j,0)
        enddo
@@ -68,9 +191,10 @@
       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,(-gcart(j,i),j=1,3)
+             i,(-gxcart(j,i),j=1,3)
           do j=1,3
             ind=ind+1
             zapas(ind)=-gxcart(j,i)
@@ -91,14 +215,17 @@
         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)
           enddo
         endif
       enddo
-      
+
+#endif      
       if(lmuca) then
        imtime=imtime+1
        if(mucadyn.gt.0) call muca_update(potE)       
@@ -152,6 +279,7 @@
       include 'mpif.h'
       integer :: ierror
       real(kind=8) :: time00
+      real(kind=8) ,allocatable, dimension(:)  :: DDM,DDU1,DDU2
 #endif
 !      include 'COMMON.SETUP'
 !      include 'COMMON.VAR'
@@ -168,10 +296,14 @@
 !#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
+#endif
       real(kind=8) :: dtdi
       real(kind=8),dimension(2*nres) :: massvec,sqreig	!(maxres2) maxres2=2*maxres
+      real(kind=8) :: relfeh,eps1,eps2
 !el      real(kind=8),dimension(:),allocatable :: Ghalf
 !el      real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf	!(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2))
 !el      real(kind=8),dimension(2*nres,2*nres) :: Gcopy	!(maxres2,maxres2)
@@ -180,11 +312,18 @@
       integer,dimension(6*nres) :: iwork	!(maxres6) maxres6=6*maxres
 !el      common /przechowalnia/ Gcopy,Ghalf
       real(kind=8) :: coeff
-      integer :: i,j,ind,ind1,k,ii,jj,m,m1,ii1,iti,nres2,ierr
+      integer :: i,j,ind,ind1,k,l,ii,jj,m,m1,ii1,iti,nres2,ierr,nind,mark
+      real(kind=8) :: ip4
+      integer :: iz,mnum
+      relfeh=1.0d-14
       nres2=2*nres
 
+      write (iout,*) "before FIVEDIAG"
+#ifndef FIVEDIAG
+      write (iout,*) "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)
+#endif
 !
 ! Set up the matrix of the (dC,dX)-->(C,X) transformation (A), the
 ! inertia matrix (Gmat) and the inverse of the inertia matrix (Ginv)
@@ -194,6 +333,144 @@
       dimen=(nct-nnt+1)+nside
       dimen1=(nct-nnt)+(nct-nnt+1)
       dimen3=dimen*3
+      write (iout,*) "nnt",nnt," nct",nct," nside",nside
+#ifdef FIVEDIAG
+       ip4=ip/4
+       DM(1)=mp/4+ip4
+       if (iabs(itype(nnt)).eq.10) then
+         DM(1)=DM(1)+msc(10)
+         ind=2
+         nind=1
+       else
+         DM(1)=DM(1)+isc(iabs(itype(nnt)))
+         DM(2)=msc(iabs(itype(nnt)))+isc(iabs(itype(nnt)))
+         ind=3
+         nind=2
+       endif	  
+       do i=nnt+1,nct-1
+!         if (iabs(itype(i,1)).eq.ntyp1) cycle
+         DM(ind)=2*ip4+mp/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)
+           ind=ind+1
+         else
+           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
+       if (nct.gt.nnt) then
+       DM(ind)=ip4+mp/4
+       if (iabs(itype(nct)).eq.10) then
+         DM(ind)=DM(ind)+msc(10)
+         nind=ind
+       else
+         DM(ind)=DM(ind)+isc(iabs(itype(nct)))
+         DM(ind+1)=msc(iabs(itype(nct)))+isc(iabs(itype(nct)))
+         nind=ind+1
+	endif
+        endif
+	
+	
+        ind=1
+        do i=nnt,nct
+	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
+	  DU1(ind)=mp/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/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
+       DMorig=DM
+       DU1orig=DU1
+       DU2orig=DU2
+       write (iout,*) "nind",nind," dimen",dimen
+       if (nind.ne.dimen) then
+         write (iout,*) "ERROR, dimensions differ"
+#ifdef MPI
+         call MPI_Finalize(ierr)
+#endif
+         stop
+       endif
+     write (iout,*) "The upper part of the five-diagonal inertia matrix"
+       do i=1,dimen
+         if (i.lt.dimen-1) then
+           write (iout,'(i3,3f10.5)') i,DM(i),DU1(i),DU2(i)
+         else if (i.eq.dimen-1) then  
+           write (iout,'(i3,3f10.5)') i,DM(i),DU1(i)
+         else
+           write (iout,'(i3,3f10.5)') i,DM(i)
+         endif 
+       enddo
+
+        call FDISYP (dimen, DM, DU1, DU2, MARK)
+
+        if (mark.eq.-1) then
+     write (iout,*) "ERROR: the inertia matrix is not positive definite"
+#ifdef MPI
+         call MPI_Finalize(ierr)
+#endif
+         stop
+        else if (mark.eq.0) then
+     write (iout,*) "ERROR: the inertia matrix is singular"
+#ifdef MPI
+         call MPI_Finalize(ierr)
+#endif
+        else if (mark.eq.1) then
+          write (iout,*) "The transformed inertia matrix"
+          do i=1,dimen
+            if (i.lt.dimen-1) then
+              write (iout,'(i3,3f10.5)') i,DM(i),DU1(i),DU2(i)
+            else if (i.eq.dimen-1) then  
+              write (iout,'(i3,3f10.5)') i,DM(i),DU1(i)
+            else
+              write (iout,'(i3,3f10.5)') i,DM(i)
+            endif 
+          enddo
+        endif
+! Diagonalization of the pentadiagonal matrix
+        allocate(DDU1(2*nres))
+        allocate(DDU2(2*nres))
+        allocate(DDM(2*nres))
+        do i=1,dimen-1
+          DDU1(i+1)=DU1orig(i)
+        enddo
+        do i=1,dimen-2
+          DDU2(i+2)=DU2orig(i)
+        enddo
+        DDM=DMorig
+        call quindibisect2(DDM,DDU1,DDU2,dimen,1,dimen,eps1,relfeh,eps2,iz,geigen)
+        if (lprn) then
+          write (iout,*) &
+           "Eigenvalues of the five-diagonal inertia matrix"
+          do i=1,dimen
+            write (iout,'(i5,f10.5)') i,geigen(i)
+          enddo
+        endif
+        deallocate(DDU1)
+        deallocate(DDU2)
+        deallocate(DDM)
+#else 
+! Old Gmatrix
 #ifdef MPI
       if (nfgtasks.gt.1) then
       time00=MPI_Wtime()
@@ -230,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
       
@@ -251,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
@@ -314,6 +597,49 @@
         write (iout,'(//a)') "Ginv"
         call matout(dimen,dimen,nres2,nres2,Ginv)
       endif
+
+      Bmat=0.0d0
+      Bmat(1,1)=1.0d0
+      j=2
+      do i=nnt,nct
+        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
+        else
+          if (i.lt.nct) then
+            Bmat(i-nnt+2,j-1)=-1
+            Bmat(i-nnt+2,j+1)=1
+            Bmat(i-nnt+1+(nct-nnt)+1,j-1)=-1
+            Bmat(i-nnt+1+(nct-nnt)+1,j)=1
+            j=j+2
+          else
+            Bmat(i-nnt+1+(nct-nnt)+1,j-1)=-1
+            Bmat(i-nnt+1+(nct-nnt)+1,j)=1
+          endif
+        endif
+      enddo
+      write (iout,*) "j",j," dimen",dimen 
+      if (lprn) then
+        write (iout,'(//a)') "Bmat"
+        call matout(dimen,dimen,nres2,nres2,Bmat)
+      endif
+      Gcopy=0.0d0
+      do i=1,dimen
+        do j=1,dimen
+            do k=1,dimen
+               do l=1,dimen
+                 Gcopy(i,j)=Gcopy(i,j)+Bmat(k,i)*Gmat(k,l)*Bmat(l,j)
+               enddo
+            enddo
+        enddo
+      enddo
+      if (lprn) then
+        write (iout,'(//a)') "Gmat-transformed"
+        call matout(dimen,dimen,nres2,nres2,Gcopy)
+      endif
 #ifdef MPI
       if (nfgtasks.gt.1) then
         myginv_ng_count=nres2*my_ng_count
@@ -388,6 +714,7 @@
         enddo
       endif
       deallocate(Gcopy)
+#endif
       return
       end subroutine setup_MD_matrices
 !-----------------------------------------------------------------------------