X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2FMD.f90;h=fa7131db20cd082bb95ba4ac4533a1952d8c1376;hb=a01e4f4abfd7680ce30209afd6cab5d52532a55f;hp=c509ee1c7d835cf8e8a402021e24702c8e9de7d0;hpb=299e2c41124d3fa8adba7244716515a2cc160ed1;p=unres4.git

diff --git a/source/unres/MD.f90 b/source/unres/MD.f90
index c509ee1..fa7131d 100644
--- a/source/unres/MD.f90
+++ b/source/unres/MD.f90
@@ -2570,6 +2570,7 @@
 !      implicit real*8 (a-h,o-z)
       use energy_data
       use random, only:anorm_distr
+      use MD_data
 !      include 'DIMENSIONS'
 !      include 'COMMON.CONTROL'
 !      include 'COMMON.VAR'
@@ -2588,7 +2589,15 @@
 !      include 'COMMON.NAMES'
 !      include 'COMMON.TIME1'
       real(kind=8) :: xv,sigv,lowb,highb  ,Ek1
-      integer :: i,j,ii,k,ind
+#ifdef FIVEDIAG
+       real(kind=8) ,allocatable, dimension(:)  :: DDU1,DDU2,DL2,DL1,xsolv,DML,rs
+       real(kind=8) :: sumx
+#ifdef DEBUG
+       real(kind=8) ,allocatable, dimension(:)  :: rsold
+       real (kind=8),allocatable,dimension(:,:) :: matold
+#endif
+#endif
+      integer :: i,j,ii,k,ind,mark,imark
 ! Generate random velocities from Gaussian distribution of mean 0 and std of KT/m 
 ! First generate velocities in the eigenspace of the G matrix
 !      write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
@@ -2602,22 +2611,303 @@
           lowb=-5*sigv
           highb=5*sigv
           d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
-!          write (iout,*) "i",i," ii",ii," geigen",geigen(i),&
-!            " d_t_work_new",d_t_work_new(ii)
+#ifdef DEBUG
+          write (iout,*) "i",i," ii",ii," geigen",geigen(i),&
+            " d_t_work_new",d_t_work_new(ii)
+#endif
         enddo
       enddo
+#ifdef DEBUG
 ! diagnostics
-!      Ek1=0.0d0
-!      ii=0
-!      do i=1,dimen
-!        do k=1,3
-!          ii=ii+1
-!          Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
-!        enddo
-!      enddo
-!      write (iout,*) "Ek from eigenvectors",Ek1
+      Ek1=0.0d0
+      ii=0
+      do i=1,dimen
+        do k=1,3
+          ii=ii+1
+          Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
+        enddo
+      enddo
+      write (iout,*) "Ek from eigenvectors",Ek1
+      write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
 ! end diagnostics
+#endif
+#ifdef FIVEDIAG
 ! Transform velocities to UNRES coordinate space
+       allocate (DL1(2*nres))
+       allocate (DDU1(2*nres))
+       allocate (DL2(2*nres))
+       allocate (DDU2(2*nres))
+       allocate (xsolv(2*nres))
+       allocate (DML(2*nres))
+       allocate (rs(2*nres))
+#ifdef DEBUG
+       allocate (rsold(2*nres))
+       allocate (matold(2*nres,2*nres))
+       matold=0
+       matold(1,1)=DMorig(1)
+       matold(1,2)=DU1orig(1)
+       matold(1,3)=DU2orig(1)
+       write (*,*) DMorig(1),DU1orig(1),DU2orig(1)
+
+      do i=2,dimen-1
+        do j=1,dimen
+          if (i.eq.j) then
+             matold(i,j)=DMorig(i)
+             matold(i,j-1)=DU1orig(i-1)
+           if (j.ge.3) then
+                 matold(i,j-2)=DU2orig(i-2)
+               endif
+
+            endif
+
+          enddo
+          do j=1,dimen-1
+            if (i.eq.j) then
+              matold(i,j+1)=DU1orig(i)
+
+             end if
+          enddo
+          do j=1,dimen-2
+            if(i.eq.j) then
+               matold(i,j+2)=DU2orig(i)
+            endif
+          enddo
+       enddo
+       matold(dimen,dimen)=DMorig(dimen)
+       matold(dimen,dimen-1)=DU1orig(dimen-1)
+       matold(dimen,dimen-2)=DU2orig(dimen-2)
+       write(iout,*) "old gmatrix"
+       call matout(dimen,dimen,2*nres,2*nres,matold)
+#endif
+      d_t_work=0.0d0
+      do i=1,dimen
+! Find the ith eigenvector of the pentadiagonal inertiq matrix
+       do imark=dimen,1,-1
+
+         do j=1,imark-1
+           DML(j)=DMorig(j)-geigen(i)
+         enddo
+         do j=imark+1,dimen
+           DML(j-1)=DMorig(j)-geigen(i)
+         enddo
+         do j=1,imark-2
+           DDU1(j)=DU1orig(j)
+         enddo
+         DDU1(imark-1)=DU2orig(imark-1)
+         do j=imark+1,dimen-1
+           DDU1(j-1)=DU1orig(j)
+         enddo
+         do j=1,imark-3
+           DDU2(j)=DU2orig(j)
+         enddo
+         DDU2(imark-2)=0.0d0
+         DDU2(imark-1)=0.0d0
+         do j=imark,dimen-3
+           DDU2(j)=DU2orig(j+1)
+         enddo 
+         do j=1,dimen-3
+           DL2(j+2)=DDU2(j)
+         enddo
+         do j=1,dimen-2
+           DL1(j+1)=DDU1(j)
+         enddo
+#ifdef DEBUG
+         write (iout,*) "DL2,DL1,DML,DDU1,DDU2"
+         write(iout,'(10f10.5)') (DL2(k),k=3,dimen-1)
+         write(iout,'(10f10.5)') (DL1(k),k=2,dimen-1)
+         write(iout,'(10f10.5)') (DML(k),k=1,dimen-1)
+         write(iout,'(10f10.5)') (DDU1(k),k=1,dimen-2)
+         write(iout,'(10f10.5)') (DDU2(k),k=1,dimen-3)
+#endif
+         rs=0.0d0
+         if (imark.gt.2) rs(imark-2)=-DU2orig(imark-2)
+         if (imark.gt.1) rs(imark-1)=-DU1orig(imark-1)
+         if (imark.lt.dimen) rs(imark)=-DU1orig(imark)
+         if (imark.lt.dimen-1) rs(imark+1)=-DU2orig(imark)
+#ifdef DEBUG
+         rsold=rs
+#endif
+!         write (iout,*) "Vector rs"
+!         do j=1,dimen-1
+!           write (iout,*) j,rs(j)
+!         enddo
+         xsolv=0.0d0
+         call FDIAG(dimen-1,DL2,DL1,DML,DDU1,DDU2,rs,xsolv,mark)
+
+         if (mark.eq.1) then
+
+#ifdef DEBUG
+! Check solution
+         do j=1,imark-1
+           sumx=-geigen(i)*xsolv(j)
+           do k=1,imark-1
+             sumx=sumx+matold(j,k)*xsolv(k)
+           enddo
+           do k=imark+1,dimen
+             sumx=sumx+matold(j,k)*xsolv(k-1)
+           enddo
+           write(iout,'(i5,3f10.5)') j,sumx,rsold(j),sumx-rsold(j)
+         enddo
+         do j=imark+1,dimen
+           sumx=-geigen(i)*xsolv(j-1)
+           do k=1,imark-1
+             sumx=sumx+matold(j,k)*xsolv(k)
+           enddo
+           do k=imark+1,dimen
+             sumx=sumx+matold(j,k)*xsolv(k-1)
+           enddo
+           write(iout,'(i5,3f10.5)') j-1,sumx,rsold(j-1),sumx-rsold(j-1)
+         enddo
+! end check
+         write (iout,*)&
+          "Solution of equations system",i," for eigenvalue",geigen(i) 
+         do j=1,dimen-1
+           write(iout,'(i5,f10.5)') j,xsolv(j) 
+         enddo
+#endif
+         do j=dimen-1,imark,-1
+           xsolv(j+1)=xsolv(j)
+         enddo
+         xsolv(imark)=1.0d0
+#ifdef DEBUG
+         write (iout,*) "Un-normalized eigenvector",i," for eigenvalue",geigen(i) 
+         do j=1,dimen
+           write(iout,'(i5,f10.5)') j,xsolv(j) 
+         enddo
+#endif
+! Normalize ith eigenvector
+         sumx=0.0d0
+         do j=1,dimen
+           sumx=sumx+xsolv(j)**2
+         enddo
+         sumx=dsqrt(sumx)
+         do j=1,dimen
+           xsolv(j)=xsolv(j)/sumx
+         enddo
+#ifdef DEBUG
+         write (iout,*) "Eigenvector",i," for eigenvalue",geigen(i) 
+         do j=1,dimen
+           write(iout,'(i5,3f10.5)') j,xsolv(j)
+         enddo
+#endif
+! All done at this point for eigenvector i, exit loop
+         exit
+
+         endif ! mark.eq.1
+
+       enddo ! imark
+
+       if (mark.ne.1) then
+         write (iout,*) "Unable to find eigenvector",i
+       endif
+
+!       write (iout,*) "i=",i
+       do k=1,3
+!         write (iout,*) "k=",k
+         do j=1,dimen
+           ind=(j-1)*3+k
+!           write(iout,*) "ind",ind," ind1",3*(i-1)+k
+           d_t_work(ind)=d_t_work(ind) &
+                            +xsolv(j)*d_t_work_new(3*(i-1)+k)
+         enddo
+       enddo
+      enddo ! i (loop over the eigenvectors)
+
+#ifdef DEBUG
+      write (iout,*) "d_t_work"
+      do i=1,3*dimen
+        write (iout,"(i5,f10.5)") i,d_t_work(i)
+      enddo
+      Ek1=0.0d0
+      ii=0
+      do i=nnt,nct
+        if (itype(i).eq.10) then
+          j=ii+3
+        else
+          j=ii+6
+        endif
+        if (i.lt.nct) then
+          do k=1,3
+!            write (iout,*) "k",k," ii+k",ii+k," ii+j+k",ii+j+k,"EK1",Ek1
+            Ek1=Ek1+0.5d0*mp*((d_t_work(ii+j+k)+d_t_work_new(ii+k))/2)**2+&
+            0.5d0*0.25d0*IP*(d_t_work(ii+j+k)-d_t_work_new(ii+k))**2
+          enddo
+        endif
+        if (itype(i).ne.10) ii=ii+3
+        write (iout,*) "i",i," itype",itype(i)," mass",msc(itype(i))
+        write (iout,*) "ii",ii
+        do k=1,3
+          ii=ii+1
+          write (iout,*) "k",k," ii",ii,"EK1",EK1
+          if (iabs(itype(i)).ne.10) Ek1=Ek1+0.5d0*Isc(iabs(itype(i)))*(d_t_work(ii)-d_t_work(ii-3))**2
+          Ek1=Ek1+0.5d0*msc(iabs(itype(i)))*d_t_work(ii)**2
+        enddo
+        write (iout,*) "i",i," ii",ii
+      enddo
+      write (iout,*) "Ek from d_t_work",Ek1
+      write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
+#endif      
+      deallocate(DDU1)
+      deallocate(DDU2)
+      deallocate(DL2)
+      deallocate(DL1)
+      deallocate(xsolv)
+      deallocate(DML)
+      deallocate(rs)
+#ifdef DEBUG
+      deallocate(matold)
+      deallocate(rsold)
+#endif
+      ind=1
+      do j=nnt,nct
+        do k=1,3
+          d_t(k,j)=d_t_work(ind)
+          ind=ind+1
+        enddo
+        if (itype(j).ne.10 .and. itype(j).ne.ntyp1) then
+          do k=1,3
+            d_t(k,j+nres)=d_t_work(ind)
+            ind=ind+1
+          enddo
+        end if
+      enddo
+#ifdef DEBUG
+      write (iout,*) "Random velocities in the Calpha,SC space"
+      do i=nnt,nct
+        write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
+      enddo
+      do i=nnt,nct
+        write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
+      enddo
+#endif
+      do j=1,3
+        d_t(j,0)=d_t(j,nnt)
+      enddo
+      do i=nnt,nct
+        if (itype(i).eq.10) then
+          do j=1,3
+            d_t(j,i)=d_t(j,i+1)-d_t(j,i)
+          enddo
+        else
+          do j=1,3
+            d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
+            d_t(j,i)=d_t(j,i+1)-d_t(j,i)
+          enddo
+        end if
+      enddo
+#ifdef DEBUG
+      write (iout,*)"Random velocities in the virtual-bond-vector space"
+      do i=nnt,nct-1
+        write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3)
+      enddo
+      do i=nnt,nct
+        write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3)
+      enddo
+      call kinetic(Ek1)
+      write (iout,*) "Ek from d_t_work",Ek1
+      write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb)
+#endif
+#else
       do k=0,2       
         do i=1,dimen
           ind=(i-1)*3+k+1
@@ -2649,11 +2939,13 @@
           enddo
         endif
       enddo
+#endif
 !      call kinetic(EK)
 !      write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",&
 !        2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
 !      call flush(iout)
       return
+#undef DEBUG
       end subroutine random_vel
 !-----------------------------------------------------------------------------
 #ifndef LANG0
@@ -5628,9 +5920,14 @@
 !      common /lagrange/
       allocate(d_t(3,0:nres2),d_a(3,0:nres2),d_t_old(3,0:nres2)) !(3,0:MAXRES2)
       allocate(d_a_work(nres6)) !(6*MAXRES)
+#ifdef FIVEDIAG
+      allocate(DM(nres2),DU1(nres2),DU2(nres2))
+      allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2))
+#else
       allocate(Gmat(nres2,nres2),A(nres2,nres2))
       if(.not.allocated(Ginv)) allocate(Ginv(nres2,nres2)) !in control: ergastulum
       allocate(Gsqrp(nres2,nres2),Gsqrm(nres2,nres2),Gvec(nres2,nres2)) !(maxres2,maxres2)
+#endif
       allocate(Geigen(nres2)) !(maxres2)
       if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
 !      common /inertia/ in io_conf: parmread