NEWGRAD added to unres Makefile and UNRES_NEWGRAD to cmake
authorCezary Czaplewski <czarek@chem.univ.gda.pl>
Fri, 14 Apr 2017 07:36:59 +0000 (09:36 +0200)
committerCezary Czaplewski <czarek@chem.univ.gda.pl>
Fri, 14 Apr 2017 07:36:59 +0000 (09:36 +0200)
-DUNRES_NEWGRAD=ON for cmake commandline, default OFF

14 files changed:
CMakeLists.txt
source/unres/CMakeLists.txt
source/unres/MD.f90
source/unres/Makefile
source/unres/REMD.f90
source/unres/cinfo.f90
source/unres/data/MD_data.f90
source/unres/energy.f90
source/unres/fdiag.f90 [new file with mode: 0644]
source/unres/fdisy.f90 [new file with mode: 0644]
source/unres/machpd.f90 [new file with mode: 0644]
source/unres/quindet2ok.F90 [new file with mode: 0644]
source/unres/quindibisectok.F90 [new file with mode: 0644]
source/unres/unres.f90

index 584d1d7..de3836d 100644 (file)
@@ -6,8 +6,8 @@ cmake_minimum_required(VERSION 2.8)
 project(UNRESPACK Fortran C)
 
 set(UNRES_MAJOR 4)
-set(UNRES_MINOR 0)
-set(UNRES_PATCH 1)
+set(UNRES_MINOR 1)
+set(UNRES_PATCH 0)
 set(UNRES_VERSION ${UNRES_MAJOR}.${UNRES_MINOR}.${UNRES_PATCH})
 
 #======================================
@@ -121,6 +121,7 @@ endif (NOT UNRES_FF)
 # Use of MPI library (default ON)
 option(UNRES_WITH_MPI "Choose whether or not to use MPI library" ON )
 
+option(UNRES_NEWGRAD "Choose whether or not to use NEWGRAD" OFF )
 
 
 #=================================
index cc0b694..2a48eda 100644 (file)
@@ -44,6 +44,13 @@ set(UNRES_MD_SRC0
         REMD.f90
 )
 
+if(UNRES_NEWGRAD)
+ set(UNRES_MD_SRC0 ${UNRES_MD_SRC0} 
+       fdisy.f90 fdiag.f90 machpd.f90 
+       quindibisectok.F90
+       quindet2ok.F90)
+endif()
+
 set(UNRES_MD_SRC1
        data/minim_data.f90 
 )
@@ -110,6 +117,11 @@ elseif(UNRES_MD_FF STREQUAL "4P")
   set(CPPFLAGS "UNRES -DISNAN -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC -DSCCORPDB" )
 endif(UNRES_MD_FF STREQUAL "GAB")
 
+if(UNRES_NEWGRAD)
+ set(CPPFLAGS "${CPPFLAGS} -DCARGRAD -DFIVEDIAG ")
+endif()
+
+
 #=========================================
 # System specific flags
 #=========================================
@@ -156,10 +168,18 @@ endif (architektura STREQUAL "64")
 #========================================
 if(UNRES_WITH_MPI) 
   # binary with mpi
-  set(UNRES_BIN "unresMD_${Fortran_COMPILER_NAME}_MPI_${UNRES_MD_FF}.exe")
+  if(UNRES_NEWGRAD)
+   set(UNRES_BIN "unresMD_${Fortran_COMPILER_NAME}_MPI_${UNRES_MD_FF}_NEWGRAD.exe")
+  else()
+   set(UNRES_BIN "unresMD_${Fortran_COMPILER_NAME}_MPI_${UNRES_MD_FF}.exe")
+  endif()
 else(UNRES_WITH_MPI)
   # binary without mpi
-  set(UNRES_BIN "unresMD_${Fortran_COMPILER_NAME}_single_${UNRES_MD_FF}.exe")
+  if(UNRES_NEWGRAD)
+   set(UNRES_BIN "unresMD_${Fortran_COMPILER_NAME}_single_${UNRES_MD_FF}_NEWGRAD.exe")
+  else()
+   set(UNRES_BIN "unresMD_${Fortran_COMPILER_NAME}_single_${UNRES_MD_FF}.exe")
+  endif()
 endif(UNRES_WITH_MPI)  
 
 #=========================================
index c509ee1..fa7131d 100644 (file)
 !      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'
 !      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
           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
           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
 !      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
index c1e9c7a..86c700c 100644 (file)
@@ -46,9 +46,9 @@ data = names.o io_units.o calc_data.o compare_data.o control_data.o \
 
 objects = xdrf/*.o \
        prng_32.o math.o random.o geometry.o md_calc.o io_base.o energy.o check_bond.o muca_md.o\
-       control.o io_config.o MPI.o minim.o \
+       control.o io_config.o MPI.o minim.o fdisy.o fdiag.o machpd.o \
        regularize.o compare.o map.o REMD.o MCM_MD.o io.o \
-       MD.o MREMD.o CSA.o unres.o
+       MD.o MREMD.o CSA.o unres.o quindibisectok.o quindet2ok.o
 
 
 #${EXE_FILE}: ${objects}
@@ -99,6 +99,15 @@ E0LL2Y: ${data} ${objects}
        ${FC} ${FFLAGS} cinfo.f90
        ${FC} ${OPT} ${data} ${objects} cinfo.o -o ${EXE_FILE}
 
+NEWGRAD: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DAMD64 -DUNRES -DISNAN -DMP -DMPI \
+        -DSPLITELE -DLANG0 -DCARGRAD -DFIVEDIAG
+NEWGRAD: EXE_FILE = ../../bin/unres_E0LL2Y_F90_EL-NEWG.exe
+NEWGRAD: ${data} ${objects}
+       cc -o compinfo compinfo.c
+       ./compinfo | true
+       ${FC} ${FFLAGS} cinfo.f90
+       ${FC} ${OPT} ${data} ${objects} cinfo.o -o ${EXE_FILE}
+
 
 xdrf/*.o:
        cd xdrf && make
@@ -221,3 +230,18 @@ CSA.o: CSA.f90
 
 unres.o: unres.f90
        ${FC} ${FFLAGS} ${CPPFLAGS} unres.f90
+
+quindibisectok.o: quindibisectok.F90
+       ${FC} ${FFLAGS} ${CPPFLAGS} quindibisectok.F90
+
+quindet2ok.o: quindet2ok.F90
+       ${FC} ${FFLAGS} ${CPPFLAGS} quindet2ok.F90
+
+
+
+
+
+
+
+
+
index edbcc8e..1deb000 100644 (file)
 !       include 'COMMON.CONTROL'
 !       include 'COMMON.MUCA'
 !      include 'COMMON.TIME1'
-       
-       integer :: i,j,ind,itime
+       integer ::i,j,ind,itime,k
        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.
 !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).ne.10 .and. itype(i).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
+           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
+         if (itype(i).eq.10 .or. itype(i).eq.ntyp1)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
+         if (itype(i).eq.10 .or. itype(i).eq.ntyp1) 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
       do i=nnt,nct
         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) 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)
           enddo
         endif
       enddo
-      
+
+#endif      
       if(lmuca) then
        imtime=imtime+1
        if(mucadyn.gt.0) call muca_update(potE)       
       include 'mpif.h'
       integer :: ierror
       real(kind=8) :: time00
+      real(kind=8) ,allocatable, dimension(:)  :: DDM,DDU1,DDU2
 #endif
 !      include 'COMMON.SETUP'
 !      include 'COMMON.VAR'
 !#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)
       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
+      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)
       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)).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)
+           ind=ind+1
+         else
+           DM(ind)=DM(ind)+isc(iabs(itype(i)))
+           DM(ind+1)=msc(iabs(itype(i)))+isc(iabs(itype(i)))
+           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)).ne.10 .and.iabs(itype((i))).ne.ntyp1) then
+          DU1(ind)=-isc(iabs(itype(i)))
+          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
+!        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
+         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()
         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
+        if (itype(i).eq.10) 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
           enddo
         enddo
       endif
-!      deallocate(Gcopy)
+      deallocate(Gcopy)
+#endif
       return
       end subroutine setup_MD_matrices
 !-----------------------------------------------------------------------------
index d38c220..22da975 100644 (file)
@@ -1,16 +1,16 @@
 ! DO NOT EDIT THIS FILE - IT HAS BEEN GENERATED BY COMPINFO.C
-! 0 40376 67
+! 0 40376 71
       subroutine cinfo
 !      include 'COMMON.IOUNITS'
       use io_units
       write(iout,*)'++++ Compile info ++++'
-      write(iout,*)'Version 0.40376 build 67'
-      write(iout,*)'compiled Fri Mar 10 14:56:02 2017'
-      write(iout,*)'compiled by emilial@piasek4'
+      write(iout,*)'Version 0.40376 build 71'
+      write(iout,*)'compiled Thu Apr 13 11:26:27 2017'
+      write(iout,*)'compiled by czarek@piasek4'
       write(iout,*)'OS name:    Linux '
-      write(iout,*)'OS release: 3.2.0-111-generic '
+      write(iout,*)'OS release: 3.2.0-124-generic '
       write(iout,*)'OS version:',&
-       ' #153-Ubuntu SMP Wed Sep 21 21:23:31 UTC 2016 '
+       ' #167-Ubuntu SMP Fri Mar 3 15:25:36 UTC 2017 '
       write(iout,*)'flags:'
       write(iout,*)'INSTALL_DIR = /users/software/mpich2-1.4.1p1_in...'
       write(iout,*)'FC= ${INSTALL_DIR}/bin/mpif90'
@@ -34,6 +34,8 @@
       write(iout,*)'4P: EXE_FILE = ../../bin/unres_4P_F90_EL.exe'
       write(iout,*)'E0LL2Y: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DAMD...'
       write(iout,*)'E0LL2Y: EXE_FILE = ../../bin/unres_E0LL2Y_F90_E...'
+      write(iout,*)'NEWGRAD: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DAM...'
+      write(iout,*)'NEWGRAD: EXE_FILE = ../../bin/unres_E0LL2Y_F90_...'
       write(iout,*)'++++ End of compile info ++++'
       return
       end subroutine cinfo
index bee0a24..9ec0a10 100644 (file)
@@ -52,6 +52,8 @@
       real(kind=8),dimension(:),allocatable :: d_a_work !(6*MAXRES)
       real(kind=8),dimension(:,:),allocatable :: Gmat,Ginv,A,&
        Gsqrp,Gsqrm,Gvec !(maxres2,maxres2)
+      real(kind=8),dimension(:),allocatable :: DM,DU1,DU2
+      real(kind=8),dimension(:),allocatable :: DMorig,DU1orig,DU2orig
       real(kind=8),dimension(:),allocatable :: Geigen !(maxres2)
       real(kind=8),dimension(:),allocatable ::vtot !(maxres2)
       logical :: reset_moment,reset_vel,rattle,RESPA
index fdf4576..da1a2f4 100644 (file)
       enddo
       return
       end subroutine check_ecart
+#ifdef CARGRAD
+!-----------------------------------------------------------------------------
+      subroutine check_ecartint
+! Check the gradient of the energy in Cartesian coordinates. 
+      use io_base, only: intout
+!      implicit real*8 (a-h,o-z)
+!      include 'DIMENSIONS'
+!      include 'COMMON.CONTROL'
+!      include 'COMMON.CHAIN'
+!      include 'COMMON.DERIV'
+!      include 'COMMON.IOUNITS'
+!      include 'COMMON.VAR'
+!      include 'COMMON.CONTACTS'
+!      include 'COMMON.MD'
+!      include 'COMMON.LOCAL'
+!      include 'COMMON.SPLITELE'
+      use comm_srutu
+!el      integer :: icall
+!el      common /srutu/ icall
+      real(kind=8),dimension(6) :: ggg,ggg1
+      real(kind=8),dimension(3) :: cc,xx,ddc,ddx,ddc1,ddcn
+      real(kind=8),dimension(6*nres) :: x,g !(maxvar) (maxvar=6*maxres)
+      real(kind=8),dimension(3) :: dcnorm_safe1,dcnorm_safe2,dxnorm_safe
+      real(kind=8),dimension(6,0:nres) :: grad_s,grad_s1 !(6,0:maxres)
+      real(kind=8),dimension(nres) :: phi_temp,theta_temp,alph_temp,omeg_temp !(maxres)
+      real(kind=8),dimension(0:n_ene) :: energia,energia1
+      integer :: uiparm(1)
+      real(kind=8) :: urparm(1)
+!EL      external fdum
+      integer :: i,j,k,nf
+      real(kind=8) :: rlambd,aincr,etot,etot1,etot11,etot12,etot2,&
+                   etot21,etot22
+      r_cut=2.0d0
+      rlambd=0.3d0
+      icg=1
+      nf=0
+      nfl=0
+      call intout
+!      call intcartderiv
+!      call checkintcartgrad
+      call zerograd
+      aincr=1.0D-5
+      write(iout,*) 'Calling CHECK_ECARTINT.'
+      nf=0
+      icall=0
+      write (iout,*) "Before geom_to_var"
+      call geom_to_var(nvar,x)
+      write (iout,*) "after geom_to_var"
+      write (iout,*) "split_ene ",split_ene
+      call flush(iout)
+      if (.not.split_ene) then
+        write(iout,*) 'Calling CHECK_ECARTINT if'
+        call etotal(energia)
+!elwrite(iout,*) 'Calling CHECK_ECARTINT if'
+        etot=energia(0)
+        write (iout,*) "etot",etot
+        call flush(iout)
+!el        call enerprint(energia)
+!elwrite(iout,*) 'Calling CHECK_ECARTINT if'
+        call flush(iout)
+        write (iout,*) "enter cartgrad"
+        call flush(iout)
+        call cartgrad
+!elwrite(iout,*) 'Calling CHECK_ECARTINT if'
+        write (iout,*) "exit cartgrad"
+        call flush(iout)
+        icall =1
+        do i=1,nres
+          write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3)
+        enddo
+        do j=1,3
+          grad_s(j,0)=gcart(j,0)
+        enddo
+!elwrite(iout,*) 'Calling CHECK_ECARTINT if'
+        do i=1,nres
+          do j=1,3
+            grad_s(j,i)=gcart(j,i)
+            grad_s(j+3,i)=gxcart(j,i)
+          enddo
+        enddo
+      else
+write(iout,*) 'Calling CHECK_ECARTIN else.'
+!- split gradient check
+        call zerograd
+        call etotal_long(energia)
+!el        call enerprint(energia)
+        call flush(iout)
+        write (iout,*) "enter cartgrad"
+        call flush(iout)
+        call cartgrad
+        write (iout,*) "exit cartgrad"
+        call flush(iout)
+        icall =1
+        write (iout,*) "longrange grad"
+        do i=1,nres
+          write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),&
+          (gxcart(j,i),j=1,3)
+        enddo
+        do j=1,3
+          grad_s(j,0)=gcart(j,0)
+        enddo
+        do i=1,nres
+          do j=1,3
+            grad_s(j,i)=gcart(j,i)
+            grad_s(j+3,i)=gxcart(j,i)
+          enddo
+        enddo
+        call zerograd
+        call etotal_short(energia)
+!el        call enerprint(energia)
+        call flush(iout)
+        write (iout,*) "enter cartgrad"
+        call flush(iout)
+        call cartgrad
+        write (iout,*) "exit cartgrad"
+        call flush(iout)
+        icall =1
+        write (iout,*) "shortrange grad"
+        do i=1,nres
+          write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),&
+          (gxcart(j,i),j=1,3)
+        enddo
+        do j=1,3
+          grad_s1(j,0)=gcart(j,0)
+        enddo
+        do i=1,nres
+          do j=1,3
+            grad_s1(j,i)=gcart(j,i)
+            grad_s1(j+3,i)=gxcart(j,i)
+          enddo
+        enddo
+      endif
+      write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
+!      do i=1,nres
+      do i=nnt,nct
+        do j=1,3
+          if (nnt.gt.1 .and. i.eq.nnt) ddc1(j)=c(j,1)
+          if (nct.lt.nres .and. i.eq.nct) ddcn(j)=c(j,nres)
+         ddc(j)=c(j,i) 
+         ddx(j)=c(j,i+nres) 
+          dcnorm_safe1(j)=dc_norm(j,i-1)
+          dcnorm_safe2(j)=dc_norm(j,i)
+          dxnorm_safe(j)=dc_norm(j,i+nres)
+        enddo
+       do j=1,3
+         c(j,i)=ddc(j)+aincr
+          if (nnt.gt.1 .and. i.eq.nnt) c(j,1)=c(j,1)+aincr
+          if (nct.lt.nres .and. i.eq.nct) c(j,nres)=c(j,nres)+aincr
+          if (i.gt.1) dc(j,i-1)=c(j,i)-c(j,i-1)
+          dc(j,i)=c(j,i+1)-c(j,i)
+          dc(j,i+nres)=c(j,i+nres)-c(j,i)
+          call int_from_cart1(.false.)
+          if (.not.split_ene) then
+            call etotal(energia1)
+            etot1=energia1(0)
+            write (iout,*) "ij",i,j," etot1",etot1
+          else
+!- split gradient
+            call etotal_long(energia1)
+            etot11=energia1(0)
+            call etotal_short(energia1)
+            etot12=energia1(0)
+          endif
+!- end split gradient
+!          write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
+         c(j,i)=ddc(j)-aincr
+          if (nnt.gt.1 .and. i.eq.nnt) c(j,1)=ddc1(j)-aincr
+          if (nct.lt.nres .and. i.eq.nct) c(j,nres)=ddcn(j)-aincr
+          if (i.gt.1) dc(j,i-1)=c(j,i)-c(j,i-1)
+          dc(j,i)=c(j,i+1)-c(j,i)
+          dc(j,i+nres)=c(j,i+nres)-c(j,i)
+          call int_from_cart1(.false.)
+          if (.not.split_ene) then
+            call etotal(energia1)
+            etot2=energia1(0)
+            write (iout,*) "ij",i,j," etot2",etot2
+           ggg(j)=(etot1-etot2)/(2*aincr)
+          else
+!- split gradient
+            call etotal_long(energia1)
+            etot21=energia1(0)
+           ggg(j)=(etot11-etot21)/(2*aincr)
+            call etotal_short(energia1)
+            etot22=energia1(0)
+           ggg1(j)=(etot12-etot22)/(2*aincr)
+!- end split gradient
+!            write (iout,*) "etot21",etot21," etot22",etot22
+          endif
+!          write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
+         c(j,i)=ddc(j)
+          if (nnt.gt.1 .and. i.eq.nnt) c(j,1)=ddc1(j)
+          if (nct.lt.nres .and. i.eq.nct) c(j,nres)=ddcn(j)
+          if (i.gt.1) dc(j,i-1)=c(j,i)-c(j,i-1)
+          dc(j,i)=c(j,i+1)-c(j,i)
+          dc(j,i+nres)=c(j,i+nres)-c(j,i)
+          dc_norm(j,i-1)=dcnorm_safe1(j)
+          dc_norm(j,i)=dcnorm_safe2(j)
+          dc_norm(j,i+nres)=dxnorm_safe(j)
+        enddo
+       do j=1,3
+         c(j,i+nres)=ddx(j)+aincr
+          dc(j,i+nres)=c(j,i+nres)-c(j,i)
+          call int_from_cart1(.false.)
+          if (.not.split_ene) then
+            call etotal(energia1)
+            etot1=energia1(0)
+          else
+!- split gradient
+            call etotal_long(energia1)
+            etot11=energia1(0)
+            call etotal_short(energia1)
+            etot12=energia1(0)
+          endif
+!- end split gradient
+         c(j,i+nres)=ddx(j)-aincr
+          dc(j,i+nres)=c(j,i+nres)-c(j,i)
+          call int_from_cart1(.false.)
+          if (.not.split_ene) then
+            call etotal(energia1)
+            etot2=energia1(0)
+           ggg(j+3)=(etot1-etot2)/(2*aincr)
+          else
+!- split gradient
+            call etotal_long(energia1)
+            etot21=energia1(0)
+           ggg(j+3)=(etot11-etot21)/(2*aincr)
+            call etotal_short(energia1)
+            etot22=energia1(0)
+           ggg1(j+3)=(etot12-etot22)/(2*aincr)
+!- end split gradient
+          endif
+!          write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
+         c(j,i+nres)=ddx(j)
+          dc(j,i+nres)=c(j,i+nres)-c(j,i)
+          dc_norm(j,i+nres)=dxnorm_safe(j)
+          call int_from_cart1(.false.)
+        enddo
+       write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)') &
+         i,(ggg(k),k=1,6),(grad_s(k,i),k=1,6),(ggg(k)/grad_s(k,i),k=1,6)
+        if (split_ene) then
+          write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)') &
+         i,(ggg1(k),k=1,6),(grad_s1(k,i),k=1,6),(ggg1(k)/grad_s1(k,i),&
+         k=1,6)
+         write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)') &
+         i,(ggg(k)+ggg1(k),k=1,6),(grad_s(k,i)+grad_s1(k,i),k=1,6),&
+         ((ggg(k)+ggg1(k))/(grad_s(k,i)+grad_s1(k,i)),k=1,6)
+        endif
+      enddo
+      return
+      end subroutine check_ecartint
+#else
 !-----------------------------------------------------------------------------
       subroutine check_ecartint
 ! Check the gradient of the energy in Cartesian coordinates. 
       enddo
       return
       end subroutine check_ecartint
+#endif
 !-----------------------------------------------------------------------------
       subroutine check_eint
 ! Check the gradient of energy in internal coordinates.
             (gxcart(j,i),j=1,3)
       enddo
 #endif
+#ifdef CARGRAD
+#ifdef DEBUG
+      write (iout,*) "CARGRAD"
+#endif
+      do i=nres,1,-1
+        do j=1,3
+          gcart(j,i)=-gcart(j,i)+gcart(j,i-1)-gxcart(j,i)
+!          gcart_new(j,i)=-gcart(j,i)+gcart(j,i-1)-gxcart(j,i)
+        enddo
+!        write (iout,'(i5,3f10.5,5x,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3), &
+!            (gcart_new(j,i),j=1,3),(gxcart(j,i),j=1,3)
+      enddo    
+! Correction: dummy residues
+        if (nnt.gt.1) then
+          do j=1,3
+!            gcart_new(j,nnt)=gcart_new(j,nnt)+gcart_new(j,1)
+            gcart(j,nnt)=gcart(j,nnt)+gcart(j,1)
+          enddo
+        endif
+        if (nct.lt.nres) then
+          do j=1,3
+!            gcart_new(j,nct)=gcart_new(j,nct)+gcart_new(j,nres)
+            gcart(j,nct)=gcart(j,nct)+gcart(j,nres)
+          enddo
+        endif
+#endif
 #ifdef TIMING
       time_cartgrad=time_cartgrad+MPI_Wtime()-time00
 #endif
diff --git a/source/unres/fdiag.f90 b/source/unres/fdiag.f90
new file mode 100644 (file)
index 0000000..26b245f
--- /dev/null
@@ -0,0 +1,296 @@
+![            {Systems with Five--Diagonal Matrices}*)                  \r
+      SUBROUTINE FDIAG (N, DL2, DL1, DM, DU1, DU2, RS, X, MARK) \r
+!                                                                       \r
+!*****************************************************************      \r
+!                                                                *      \r
+!     Solving a system of linear equations                       *      \r
+!                      A * X = RS                                *      \r
+!     with a five-diagonal, strongly nonsingular matrix A via    *      \r
+!     Gauss algorithm without pivoting.                          *      \r
+!     The matrix A is given as five N-vectors DL2, DL1, DM, DU1  *      \r
+!     and DU2. The linear system has the form:                   *      \r
+!                                                                *      \r
+!     DM(1)*X(1)+DU1(1)*X(2)+DU2(1)*X(3)             = RS(1)     *      \r
+!     DL1(2)*X(1)+DM(2)*X(2)+DU1(2)*X(3)+DU2(2)*X(4) = RS(2)     *      \r
+!                                                                *      \r
+!     DL2(I)*X(I-2)+DL1(I)*X(I-1)+                               *      \r
+!           +DM(I)*X(I)+DU1(I)*X(I+1)+DU2(I)*X(I+2)  = RS(I)     *      \r
+!            for I = 3, ..., N - 2, and                          *      \r
+!                                                                *      \r
+!     DL2(N-1)*X(N-3)+DL1(N-1)*X(N-2)+                           *      \r
+!             +DM(N-1)*X(N-1)+DU1(N-1)+X(N)          = RS(N-1)   *      \r
+!     DL2(N)*X(N-2)+DL1(N)*X(N-1)+DM(N)*X(N)         = RS(N)     *      \r
+!                                                                *      \r
+!                                                                *      \r
+!                                                                *      \r
+!     INPUT PARAMETERS:                                          *      \r
+!     =================                                          *      \r
+!     N     : number of equations; N > 3                         *      \r
+!     DL2   : N-vector DL2(1:N); second lower co-diagonal        *      \r
+!             DL2(3), DL2(4), ... , DL2(N)                       *      \r
+!     DL1   : N-vector DL1(1:N); lower co-diagonal               *      \r
+!             DL1(2), DL1(3), ... , DL1(N)                       *      \r
+!     DM    : N-vector DM(1:N); main diagonal                    *      \r
+!             DM(1), DM(2), ... , DM(N)                          *      \r
+!     DU1   : N-vector DU1(1:N); upper co-diagonal               *      \r
+!             DU1(1), DU1(2), ... , DU1(N-1)                     *      \r
+!     DU2   : N-vector DU2(1:N); second upper co-diagonal        *      \r
+!             DU2(1), DU2(2), ... , DU2(N-2)                     *      \r
+!     RS    : N-vector RS(1:N); the right hand side of the       *      \r
+!             linear system                                      *      \r
+!                                                                *      \r
+!                                                                *      \r
+!     OUTPUT PARAMETERS:                                         *      \r
+!     ==================                                         *      \r
+!     DL2   :) overwritten with auxiliary vectors defining the   *      \r
+!     DL1   :) factorization of the cyclically tridiagonal       *      \r
+!     DM    :) matrix A                                          *      \r
+!     DU1   :)                                                   *      \r
+!     DU2   :)                                                   *      \r
+!     X     : N-vector X(1:N); containing the solution of the    *      \r
+!             the system of equations                            *      \r
+!     MARK  : error parameter                                    *      \r
+!             MARK=-1 : condition N > 3 is not satisfied         *      \r
+!             MARK= 0 : numerically the matrix A is not strongly *      \r
+!                       nonsingular                              *      \r
+!             MARK= 1 : everything is o.k.                       *      \r
+!                                                                *      \r
+!     NOTE: if MARK = 1, the determinant of A is given by:       *      \r
+!                DET A = DM(1) * DM(2) * ... * DM(N)             *      \r
+!                                                                *      \r
+!----------------------------------------------------------------*      \r
+!                                                                *      \r
+!  subroutines required: FDIAGP, FDIAGS, MACHPD                  *      \r
+!                                                                *      \r
+!*****************************************************************      \r
+!                                                                *      \r
+!  author   : Gisela Engeln-Muellges                             *      \r
+!  date     : 05.06.1988                                         *      \r
+!  source   : FORTRAN 77                                         *      \r
+!                                                                *      \r
+!*****************************************************************      \r
+!                                                                       \r
+      IMPLICIT DOUBLEPRECISION (A - H, O - Z) \r
+      DOUBLEPRECISION DL1 (1:N), DL2 (1:N), DM (1:N) \r
+      DOUBLEPRECISION DU1 (1:N), DU2 (1:N), RS (1:N), X (1:N) \r
+      MARK = - 1 \r
+      IF (N.LT.4) RETURN \r
+!                                                                       \r
+!  Factor the matrix A                                                  \r
+!                                                                       \r
+      CALL FDIAGP (N, DL2, DL1, DM, DU1, DU2, MARK) \r
+!                                                                       \r
+!  if MARK = 1, update and bachsubstitute                               \r
+!                                                                       \r
+      IF (MARK.EQ.1) THEN \r
+         CALL FDIAGS (N, DL2, DL1, DM, DU1, DU2, RS, X) \r
+      ENDIF \r
+      RETURN \r
+      END SUBROUTINE FDIAG                          \r
+!                                                                       \r
+!                                                                       \r
+      SUBROUTINE FDIAGP (N, DL2, DL1, DM, DU1, DU2, MARK) \r
+!                                                                       \r
+!*****************************************************************      \r
+!                                                                *      \r
+!     Factor a five-diagonal, strongly nonsingular matrix A      *      \r
+!     that is defined by the five N-vectors DL2, DL1, DM, DU1    *      \r
+!     and DU2, into its triangular factors  L * R  by applying   *      \r
+!     Gaussian elimination specialized for five-diagonal matrices*      \r
+!     (without pivoting).                                        *      \r
+!                                                                *      \r
+!                                                                *      \r
+!     INPUT PARAMETERS:                                          *      \r
+!     =================                                          *      \r
+!     N     : number of equations; N > 3                         *      \r
+!     DL2   : N-vector DL2(1:N); second lower co-diagonal        *      \r
+!             DL2(3), DL2(4), ... , DL2(N)                       *      \r
+!     DL1   : N-vector DL1(1:N); lower co-diagonal               *      \r
+!             DL1(2), DL1(3), ... , DL1(N)                       *      \r
+!     DM    : N-vector DM(1:N); main diagonal                    *      \r
+!             DM(1), DM(2), ... , DM(N)                          *      \r
+!     DU1   : N-vector DU1(1:N); upper co-diagonal               *      \r
+!             DU1(1), DU1(2), ... , DU1(N-1)                     *      \r
+!     DU2   : N-vector DU2(1:N); second upper co-diagonal        *      \r
+!             DU2(1), DU2(2), ... , DU2(N-2)                     *      \r
+!                                                                *      \r
+!                                                                *      \r
+!     OUTPUT PARAMETERS:                                         *      \r
+!     ==================                                         *      \r
+!     DL2   :) overwritten with auxiliary vectors that define    *      \r
+!     DL1   :) the factors of the five-diagonal matrix A;        *      \r
+!     DM    :) the three co-diagonals of the lower triangular    *      \r
+!     DU1   :) matrix L are stored in the vectors DL2, DL1 and   *      \r
+!     DU2   :) DM. The two co-diagonals of the unit upper        *      \r
+!              triangular matrix R are stored in the vectors DU1 *      \r
+!              and DU2, its diagonal elements each have the      *      \r
+!              value  1.                                         *      \r
+!     MARK  : error parameter                                    *      \r
+!             MARK=-1 : condition N > 3 is violated              *      \r
+!             MARK= 0 : numerically the matrix is not strongly   *      \r
+!                       nonsingular                              *      \r
+!             MARK= 1 : everything is o.k.                       *      \r
+!                                                                *      \r
+!----------------------------------------------------------------*      \r
+!                                                                *      \r
+!  subroutines required: MACHPD                                  *      \r
+!                                                                *      \r
+!*****************************************************************      \r
+!                                                                *      \r
+!  author   : Gisela Engeln-Muellges                             *      \r
+!  date     : 05.06.1988                                         *      \r
+!  source   : FORTRAN 77                                         *      \r
+!                                                                *      \r
+!*****************************************************************      \r
+!                                                                       \r
+      IMPLICIT DOUBLEPRECISION (A - H, O - Z) \r
+      DOUBLEPRECISION DL2 (1:N), DL1 (1:N), DM (1:N), DU1 (1:N),        &\r
+      DU2 (1:N)                                                         \r
+!                                                                       \r
+!  testing whether N > 3                                                \r
+!                                                                       \r
+      MARK = - 1 \r
+      IF (N.LT.4) RETURN \r
+!                                                                       \r
+!  calculating the machine constant                                     \r
+!                                                                       \r
+      FMACHP = 1.0D0 \r
+   10 FMACHP = 0.5D0 * FMACHP \r
+      IF (MACHPD (1.0D0 + FMACHP) .EQ.1) GOTO 10 \r
+      FMACHP = FMACHP * 2.0D0 \r
+!                                                                       \r
+!  determining relative error bounds                                    \r
+!                                                                       \r
+      EPS = 4.0D0 * FMACHP \r
+!                                                                       \r
+!  initializing the undefined vector components                         \r
+!                                                                       \r
+      DL2 (1) = 0.0D0 \r
+      DL2 (2) = 0.0D0 \r
+      DL1 (1) = 0.0D0 \r
+      DU1 (N) = 0.0D0 \r
+      DU2 (N - 1) = 0.0D0 \r
+      DU2 (N) = 0.0D0 \r
+!                                                                       \r
+!  factoring the matrix A while checking for strong nonsingularity      \r
+!  for N=1, 2                                                           \r
+!                                                                       \r
+      ROW = DABS (DM (1) ) + DABS (DU1 (1) ) + DABS (DU2 (1) ) \r
+      IF (ROW.EQ.0.0D0) THEN \r
+         MARK = 0 \r
+         RETURN \r
+      ENDIF \r
+      D = 1.0D0 / ROW \r
+      IF (DABS (DM (1) ) * D.LE.EPS) THEN \r
+         MARK = 0 \r
+         RETURN \r
+      ENDIF \r
+      DU1 (1) = DU1 (1) / DM (1) \r
+      DU2 (1) = DU2 (1) / DM (1) \r
+      ROW = DABS (DL1 (2) ) + DABS (DM (2) ) + DABS (DU1 (2) ) + DABS ( &\r
+      DU2 (2) )                                                         \r
+      IF (ROW.EQ.0.0D0) THEN \r
+         MARK = 0 \r
+         RETURN \r
+      ENDIF \r
+      D = 1.0D0 / ROW \r
+      DM (2) = DM (2) - DL1 (2) * DU1 (1) \r
+      IF (DABS (DM (2) ) * D.LE.EPS) THEN \r
+         MARK = 0 \r
+         RETURN \r
+      ENDIF \r
+      DU1 (2) = (DU1 (2) - DL1 (2) * DU2 (1) ) / DM (2) \r
+      DU2 (2) = DU2 (2) / DM (2) \r
+!                                                                       \r
+!  factoring A while checking for strong nonsingularity of A            \r
+!                                                                       \r
+      DO 20 I = 3, N, 1 \r
+         ROW = DABS (DL2 (I) ) + DABS (DL1 (I) ) + DABS (DM (I) )       &\r
+         + DABS (DU1 (I) ) + DABS (DU2 (I) )                            \r
+         IF (ROW.EQ.0.0D0) THEN \r
+            MARK = 0 \r
+            RETURN \r
+         ENDIF \r
+         D = 1.0D0 / ROW \r
+         DL1 (I) = DL1 (I) - DL2 (I) * DU1 (I - 2) \r
+         DM (I) = DM (I) - DL2 (I) * DU2 (I - 2) - DL1 (I) * DU1 (I - 1) \r
+         IF (DABS (DM (I) ) * D.LE.EPS) THEN \r
+            MARK = 0 \r
+            RETURN \r
+         ENDIF \r
+         IF (I.LT.N) THEN \r
+            DU1 (I) = (DU1 (I) - DL1 (I) * DU2 (I - 1) ) / DM (I) \r
+         ENDIF \r
+         IF (I.LT. (N - 1) ) THEN \r
+            DU2 (I) = DU2 (I) / DM (I) \r
+         ENDIF \r
+   20 END DO \r
+      MARK = 1 \r
+      RETURN \r
+      END SUBROUTINE FDIAGP                         \r
+!                                                                       \r
+!                                                                       \r
+      SUBROUTINE FDIAGS (N, DL2, DL1, DM, DU1, DU2, RS, X) \r
+!                                                                       \r
+!*****************************************************************      \r
+!                                                                *      \r
+!     Solving a linear system of equations                       *      \r
+!                A * X = RS                                      *      \r
+!     for a five-diagonal, strongly nonsingular matrix A, once   *      \r
+!     the factor matrices L * R have been calculated by          *      \r
+!     SUBROUTINE FDIAGP. Here they are used as input arrays and  *      \r
+!     they are stored in the five N-vectors DL2, DL1, DM, DU1    *      \r
+!     and DU2.                                                   *      \r
+!                                                                *      \r
+!                                                                *      \r
+!     INPUT PARAMETERS:                                          *      \r
+!     =================                                          *      \r
+!     N     : number of equations; N > 3                         *      \r
+!     DL2   : N-vector DL2(1:N); ) lower triangular matrix L     *      \r
+!     DL1   : N-vector DL1(1:N); ) including the diagonal        *      \r
+!     DM    : N-vector DM(1:N);  ) elements                      *      \r
+!                                                                *      \r
+!     DU1   : N-vector DU1(1:N); ) unit upper triangular matrix  *      \r
+!     DU2   : N-vector DU2(1:N); ) R without its unit diagonal   *      \r
+!                                   elements                     *      \r
+!     RS    : N-vector RS1(1:N); right side of the linear system *      \r
+!                                                                *      \r
+!                                                                *      \r
+!     OUTPUT PARAMETERS:                                         *      \r
+!     ==================                                         *      \r
+!     X     : N-vector X(1:N); the solution of the linear system *      \r
+!                                                                *      \r
+!----------------------------------------------------------------*      \r
+!                                                                *      \r
+!  subroutines required: none                                    *      \r
+!                                                                *      \r
+!*****************************************************************      \r
+!                                                                *      \r
+!  author   : Gisela Engeln-Muellges                             *      \r
+!  date     : 05.06.1988                                         *      \r
+!  source   : FORTRAN 77                                         *      \r
+!                                                                *      \r
+!*****************************************************************      \r
+!                                                                       \r
+      IMPLICIT DOUBLEPRECISION (A - H, O - Z) \r
+      DOUBLEPRECISION DL2 (1:N), DL1 (1:N), DM (1:N) \r
+      DOUBLEPRECISION DU1 (1:N), DU2 (1:N), RS (1:N), X (1:N) \r
+!                                                                       \r
+!  updating                                                             \r
+!                                                                       \r
+      RS (1) = RS (1) / DM (1) \r
+      RS (2) = (RS (2) - DL1 (2) * RS (1) ) / DM (2) \r
+      DO 10 I = 3, N \r
+         RS (I) = (RS (I) - DL2 (I) * RS (I - 2) - DL1 (I) * RS (I - 1) &\r
+         ) / DM (I)                                                     \r
+   10 END DO \r
+!                                                                       \r
+!  backsubstitution                                                     \r
+!                                                                       \r
+      X (N) = RS (N) \r
+      X (N - 1) = RS (N - 1) - DU1 (N - 1) * X (N) \r
+      DO 20 I = N - 2, 1, - 1 \r
+         X (I) = RS (I) - DU1 (I) * X (I + 1) - DU2 (I) * X (I + 2) \r
+   20 END DO \r
+      RETURN \r
+      END SUBROUTINE FDIAGS                         \r
diff --git a/source/unres/fdisy.f90 b/source/unres/fdisy.f90
new file mode 100644 (file)
index 0000000..de6f33d
--- /dev/null
@@ -0,0 +1,300 @@
+![  {Systems with Five--Diagonal Symmetric Matrices}                    \r
+![  {Systems with Five--Diagonal Symmetric Matrices}*)                  \r
+      SUBROUTINE FDISY (N, DM, DU1, DU2, RS, X, MARK) \r
+!                                                                       \r
+!*****************************************************************      \r
+!                                                                *      \r
+!  Solving a system of linear equations                          *      \r
+!                     A * X = RS                                 *      \r
+!  for a five-diagonal, symmetric and strongly nonsingular       *      \r
+!  matrix A. The matrix A is given by the three N-vectors DM,    *      \r
+!  DU1 and DU2. The system of equations has the form :           *      \r
+!                                                                *      \r
+!  DM(1)*X(1) + DU1(1)*X(2) + DU2(1)*X(3)               = RS(1)  *      \r
+!  DU1(1)*X(1) + DM(2)*X(2) + DU1(2)*X(3) + DU2(2)*X(4) = RS(2)  *      \r
+!                                                                *      \r
+!  DU2(I-2)*X(I-2) + DU1(I-1)*X(I-1) + DM(I)*X(I) +              *      \r
+!                       + DU1(I)*X(I+1) + DU2(I)*X(I+2) = RS(I)  *      \r
+!             for I = 3, ..., N - 2, and                         *      \r
+!                                                                *      \r
+!  DU2(N-3)*X(N-2) + DU1(N-2)*X(N-1) + DM(N-1)*X(N-1) +          *      \r
+!                                       + DU1(N-1)*X(N) = RS(N-1)*      \r
+!  DU2(N-2)*X(N-2) + OD(N-1)*X(N-1) + DM(N)*X(N)        = RS(N)  *      \r
+!                                                                *      \r
+!                                                                *      \r
+!                                                                *      \r
+!  INPUT PARAMETERS:                                             *      \r
+!  =================                                             *      \r
+!  N    : number of equations, N > 3                             *      \r
+!  DM   : N-vector DM(1:N); main diagonal of A                   *      \r
+!         DM(1), DM(2), ... , DM(N)                              *      \r
+!  DU1  : N-vector DU1(1:N); co-diagonal of A                    *      \r
+!         DU1(1), DU1(2), ... , DU1(N-1)                         *      \r
+!  DU2  : N-vector DU2(1:N); second co-diagonal of A             *      \r
+!         DU2(1), DU2(2), ... , DU2(N-2)                         *      \r
+!  RS   : N-vector RS(1:N); the right hand side                  *      \r
+!                                                                *      \r
+!                                                                *      \r
+!  OUTPUT PARAMETERS:                                            *      \r
+!  ==================                                            *      \r
+!  DM   :)                                                       *      \r
+!  DU1  :) overwritten with intermediate quantities              *      \r
+!  DU2  :)                                                       *      \r
+!  RS   :)                                                       *      \r
+!  X    : N-vector X(1:N) containing the solution vector         *      \r
+!  MARK : error parameter                                        *      \r
+!         MARK=-2 : condition N > 3 is not satisfied             *      \r
+!         MARK=-1 : A is strongly nonsingular, but not positive  *      \r
+!                   definite                                     *      \r
+!         MARK= 0 : numerically the matrix A is not strongly     *      \r
+!                   nonsingular                                  *      \r
+!         MARK= 1 : A is positive definite                       *      \r
+!                                                                *      \r
+!  NOTE: If MARK = +/- 1, then the determinant of A is:          *      \r
+!           DET A = DM(1) * DM(2) * ... * DM(N)                  *      \r
+!                                                                *      \r
+!----------------------------------------------------------------*      \r
+!                                                                *      \r
+!  subroutines required: FDISYP, FDISYS, MACHPD                  *      \r
+!                                                                *      \r
+!*****************************************************************      \r
+!                                                                *      \r
+!  authors  : Gisela Engeln-Muellges                             *      \r
+!  date     : 01.07.1992                                         *      \r
+!  source   : FORTRAN 77                                         *      \r
+!                                                                *      \r
+!*****************************************************************      \r
+!                                                                       \r
+      IMPLICIT DOUBLEPRECISION (A - H, O - Z) \r
+      DOUBLEPRECISION DM (1:N), DU1 (1:N), DU2 (1:N), RS (1:N), X (1:N) \r
+      MARK = - 2 \r
+      IF (N.LT.4) RETURN \r
+!                                                                       \r
+!  Factorization of the matrix A                                        \r
+!                                                                       \r
+      CALL FDISYP (N, DM, DU1, DU2, MARK) \r
+!                                                                       \r
+!  if MARK = +/- 1 , update and backsubstitute                          \r
+!                                                                       \r
+      IF (MARK.EQ.1) THEN \r
+         CALL FDISYS (N, DM, DU1, DU2, RS, X) \r
+      ENDIF \r
+      RETURN \r
+      END SUBROUTINE FDISY                          \r
+!                                                                       \r
+!                                                                       \r
+      SUBROUTINE FDISYP (N, DM, DU1, DU2, MARK) \r
+!                                                                       \r
+!*****************************************************************      \r
+!                                                                *      \r
+!  Factor a five-diagonal, symmetric and strongly nonsingular    *      \r
+!  matrix A, that is given by the three N-vectors DM, DU1 and    *      \r
+!  DU2, into its Cholesky factors A =  R(TRANSP) * D * R  by     *      \r
+!  applying the root-free Cholesky method for five-diagonal      *      \r
+!  matrices. The form of the linear system is identical with     *      \r
+!  the one in SUBROUTINE FDISY.                                  *      \r
+!                                                                *      \r
+!                                                                *      \r
+!  INPUT PARAMETERS:                                             *      \r
+!  =================                                             *      \r
+!  N    : number of equations, N > 3                             *      \r
+!  DM   : N-vector DM(1:N); main diagonal of A                   *      \r
+!         DM(1), DM(2), ... , DM(N)                              *      \r
+!  DU1  : N-vector DU1(1:N); upper co-diagonal of A              *      \r
+!         DU1(1), DU1(2), ... , DU1(N-1)                         *      \r
+!  DU2  : N-vector DU2(1:N); second upper co-diagonal of A       *      \r
+!         DU2(1), DU2(2), ... , DU2(N-2);                        *      \r
+!         due to symmetry the lower co-diagonals do not need to  *      \r
+!         be stored separately.                                  *      \r
+!                                                                *      \r
+!                                                                *      \r
+!  OUTPUT PARAMETERS:                                            *      \r
+!  ==================                                            *      \r
+!  DM   :) overwritten with auxiliary vectors containing the     *      \r
+!  DU1  :) Cholesky factors of A. The co-diagonals of the unit   *      \r
+!  DU2  :) upper tridiagonal matrix R are stored in DU1 and DU2, *      \r
+!          the diagonal matrix D in DM.                          *      \r
+!  MARK : error parameter                                        *      \r
+!         MARK=-2 : condition N > 3 is not satisfied             *      \r
+!         MARK=-1 : A is strongly nonsingular, but not positive  *      \r
+!                   definite                                     *      \r
+!         MARK= 0 : numerically the matrix is not strongly       *      \r
+!                   nonsingular                                  *      \r
+!         MARK= 1 : A is positive definite                       *      \r
+!                                                                *      \r
+!  NOTE : If MARK = +/-1, then the inertia of A, i. e., the      *      \r
+!         number of positive and negative eigenvalues of A,      *      \r
+!         is the same as the number of positive and negative     *      \r
+!         numbers among the components of DM.                    *      \r
+!                                                                *      \r
+!----------------------------------------------------------------*      \r
+!                                                                *      \r
+!  subroutines required: MACHPD                                  *      \r
+!                                                                *      \r
+!*****************************************************************      \r
+!                                                                *      \r
+!  authors  : Gisela Engeln-Muellges                             *      \r
+!  date     : 01.07.1988                                         *      \r
+!  source   : FORTRAN 77                                         *      \r
+!                                                                *      \r
+!*****************************************************************      \r
+!                                                                       \r
+      IMPLICIT DOUBLEPRECISION (A - H, O - Z) \r
+      DOUBLEPRECISION DM (1:N), DU1 (1:N), DU2 (1:N) \r
+!                                                                       \r
+!   calculating the machine constant                                    \r
+!                                                                       \r
+      FMACHP = 1.0D0 \r
+   10 FMACHP = 0.5D0 * FMACHP \r
+      IF (MACHPD (1.0D0 + FMACHP) .EQ.1) GOTO 10 \r
+      FMACHP = FMACHP * 2.0D0 \r
+!                                                                       \r
+!   determining the relative error bound                                \r
+!                                                                       \r
+      EPS = 4.0D0 * FMACHP \r
+!                                                                       \r
+!   checking for N > 3                                                  \r
+!                                                                       \r
+      MARK = - 2 \r
+      IF (N.LT.4) RETURN \r
+      DU1 (N) = 0.0D0 \r
+      DU2 (N) = 0.0D0 \r
+      DU2 (N - 1) = 0.0D0 \r
+!                                                                       \r
+!   checking for strong nonsingularity of the matrix A for N=1          \r
+!                                                                       \r
+      ROW = DABS (DM (1) ) + DABS (DU1 (1) ) + DABS (DU2 (1) ) \r
+      IF (ROW.EQ.0.0D0) THEN \r
+         MARK = 0 \r
+         RETURN \r
+      ENDIF \r
+      D = 1.0D0 / ROW \r
+      IF (DM (1) .LT.0.0D0) THEN \r
+         MARK = - 1 \r
+         RETURN \r
+      ELSEIF (DABS (DM (1) ) * D.LE.EPS) THEN \r
+         MARK = 0 \r
+         RETURN \r
+      ENDIF \r
+!                                                                       \r
+!   factoring A while checking for strong nonsingularity                \r
+!                                                                       \r
+      DUMMY = DU1 (1) \r
+      DU1 (1) = DU1 (1) / DM (1) \r
+      DUMMY1 = DU2 (1) \r
+      DU2 (1) = DU2 (1) / DM (1) \r
+      ROW = DABS (DUMMY) + DABS (DM (2) ) + DABS (DU1 (2) ) + DABS (DU2 &\r
+      (2) )                                                             \r
+      IF (ROW.EQ.0.0D0) THEN\r
+         MARK = 0\r
+         RETURN \r
+      ENDIF \r
+      D = 1.0D0 / ROW \r
+      DM (2) = DM (2) - DUMMY * DU1 (1) \r
+      IF (DM (2) .LT.0.0D0) THEN \r
+         MARK = - 1\r
+         RETURN \r
+      ELSEIF (DABS (DM (2) ) .LE.EPS) THEN \r
+         MARK = 0\r
+         RETURN \r
+      ENDIF \r
+      DUMMY = DU1 (2) \r
+      DU1 (2) = (DU1 (2) - DUMMY1 * DU1 (1) ) / DM (2) \r
+      DUMMY2 = DU2 (2) \r
+      DU2 (2) = DU2 (2) / DM (2) \r
+      DO 20 I = 3, N, 1 \r
+         ROW = DABS (DUMMY1) + DABS (DUMMY) + DABS (DM (I) ) + DABS (   &\r
+         DU1 (I) ) + DABS (DU2 (I) )                                    \r
+         IF (ROW.EQ.0.0D0) THEN \r
+            MARK = 0\r
+            RETURN \r
+         ENDIF \r
+         D = 1.0D0 / ROW \r
+         DM (I) = DM (I) - DM (I - 1) * DU1 (I - 1) * DU1 (I - 1)       &\r
+         - DUMMY1 * DU2 (I - 2)                                         \r
+         IF (DM (I) .LT.0.0D0) THEN \r
+            MARK = - 1\r
+            RETURN \r
+         ELSEIF (DABS (DM (I) ) * D.LE.EPS) THEN \r
+            MARK = 0 \r
+            RETURN \r
+         ENDIF \r
+         IF (I.LT.N) THEN \r
+            DUMMY = DU1 (I) \r
+            DU1 (I) = (DU1 (I) - DUMMY2 * DU1 (I - 1) ) / DM (I) \r
+            DUMMY1 = DUMMY2 \r
+         ENDIF \r
+         IF (I.LT.N - 1) THEN \r
+            DUMMY2 = DU2 (I) \r
+            DU2 (I) = DU2 (I) / DM (I) \r
+         ENDIF \r
+   20 END DO \r
+      MARK = 1 \r
+      RETURN \r
+      END SUBROUTINE FDISYP                         \r
+!                                                                       \r
+!                                                                       \r
+      SUBROUTINE FDISYS (N, DM, DU1, DU2, RS, X) \r
+!                                                                       \r
+!*****************************************************************      \r
+!                                                                *      \r
+!  Solving a linear system of equations                          *      \r
+!               A * X = RS                                       *      \r
+!  for a five-diagonal, symmetric and strongly nonsingular       *      \r
+!  matrix A, once its Cholesky factors have been calculated by   *      \r
+!  SUBROUTINE FDISYP. Here the factors of A are used as input    *      \r
+!  arrays and they are stored in the three N-vectors DM, DU1     *      \r
+!  and DU2.                                                      *      \r
+!                                                                *      \r
+!                                                                *      \r
+!  INPUT PARAMETER:                                              *      \r
+!  ================                                              *      \r
+!  N    : number of equations, N > 3                             *      \r
+!  DM   : N-vector DM(1:N);  diagonal matrix D                   *      \r
+!  DU1  : N-vector DM(1:N); ) co-diagonals of the upper          *      \r
+!  DU2  : N-vector DM(1:N); ) triangular  matrix R               *      \r
+!  RS   : N-vector DM(1:N); the right hand side                  *      \r
+!                                                                *      \r
+!                                                                *      \r
+!  OUTPUT PARAMETER:                                             *      \r
+!  =================                                             *      \r
+!  X    : N-vector X(1:N) containing the solution vector         *      \r
+!                                                                *      \r
+!----------------------------------------------------------------*      \r
+!                                                                *      \r
+!  subroutines required: none                                    *      \r
+!                                                                *      \r
+!*****************************************************************      \r
+!                                                                *      \r
+!  author   : Gisela Engeln-Muellges                             *      \r
+!  date     : 29.04.1988                                         *      \r
+!  source   : FORTRAN 77                                         *      \r
+!                                                                *      \r
+!*****************************************************************      \r
+!                                                                       \r
+      IMPLICIT DOUBLEPRECISION (A - H, O - Z) \r
+      DOUBLEPRECISION DM (1:N), DU1 (1:N), DU2 (1:N), RS (1:N), X (1:N) \r
+!                                                                       \r
+!  updating                                                             \r
+!                                                                       \r
+      DUMMY1 = RS (1) \r
+      RS (1) = DUMMY1 / DM (1) \r
+      DUMMY2 = RS (2) - DU1 (1) * DUMMY1 \r
+      RS (2) = DUMMY2 / DM (2) \r
+      DO 10 I = 3, N, 1 \r
+         DUMMY1 = RS (I) - DU1 (I - 1) * DUMMY2 - DU2 (I - 2) * DUMMY1 \r
+         RS (I) = DUMMY1 / DM (I) \r
+         DUMMY3 = DUMMY2 \r
+         DUMMY2 = DUMMY1 \r
+         DUMMY1 = DUMMY3 \r
+   10 END DO \r
+!                                                                       \r
+!  backsubstitution                                                     \r
+!                                                                       \r
+      X (N) = RS (N) \r
+      X (N - 1) = RS (N - 1) - DU1 (N - 1) * X (N) \r
+      DO 20 I = N - 2, 1, - 1 \r
+         X (I) = RS (I) - DU1 (I) * X (I + 1) - DU2 (I) * X (I + 2) \r
+   20 END DO \r
+      RETURN \r
+      END SUBROUTINE FDISYS                         \r
diff --git a/source/unres/machpd.f90 b/source/unres/machpd.f90
new file mode 100644 (file)
index 0000000..a297c99
--- /dev/null
@@ -0,0 +1,6 @@
+      INTEGER FUNCTION MACHPD (X) \r
+      DOUBLEPRECISION X \r
+      MACHPD = 0 \r
+      IF (1.0D0.LT.X) MACHPD = 1 \r
+      RETURN \r
+      END FUNCTION MACHPD                           \r
diff --git a/source/unres/quindet2ok.F90 b/source/unres/quindet2ok.F90
new file mode 100644 (file)
index 0000000..23e6a7b
--- /dev/null
@@ -0,0 +1,74 @@
+       subroutine quindet2(c,b,d,n,mu,no) 
+        real (kind=8) :: mu ,x
+        real (kind=8) ,dimension(n) :: c,b,d
+        real (kind=8) ,dimension(3,5) :: r
+        integer ::i,j,k,l,ll,rr,w,n,no
+! the arrays c,b i d contain the diagonal subdiagonal and subsubdiagonal  elements of a simmetric quindiagonalmatrix. the value of the output parameter   no is the number of eigenvalues greater than mu;
+        r(2,1)=c(1)-mu
+        r(2,2)=b(2)
+        r(3,1)=b(2)
+        r(2,3)=d(3)
+        r(3,2)=c(2)-mu
+        r(3,3)=b(3)
+        r(3,4)=d(4)
+        r(2,4)=0
+        r(2,5)=0
+        r(3,5)=0
+        if (r(2,1).ge.0) then
+           no=1
+        else 
+           no=0
+        endif
+! k count the major stages
+        do k=2,n 
+          rr=1
+          do i=2,1,-1
+           if (k.gt.i) then
+             w=3-i
+!interchange row w and row 3 if necessary
+             if (dabs(r(3,1)).gt. dabs(r(w,1))) then
+               do j=1,i+3 
+                 x=r(3,j)
+                 r(3,j)=r(w,j)
+                 r(w,j)=x
+               enddo
+               if (r(3,1).ge.0.eqv.r(w,1).ge.0) rr=-rr
+             endif
+!elimination of subdiagonal elements
+             if (r(w,1).eq.0) then 
+               x=0
+             else
+               x=r(3,1)/r(w,1)
+             endif
+             do j=2,i+3
+               r(3,j-1)=r(3,j)-x*r(w,j)
+             enddo
+             r(3,i+3)=0
+           endif
+          enddo ! i
+        
+          if (r(3,1).lt.0) rr=-rr
+          if (rr.gt.0) no=no+1
+! update elements of array r
+          do i = 1,2 
+            do j = 1,5
+              r(i,j)=r(i+1,j)
+            enddo
+          enddo
+          if ((k+1).le.n) then 
+            r(3,1)=d(k+1)
+            r(3,2)=b(k+1)
+            r(3,3)=c(k+1)-mu
+          endif
+          if ((k+2).le.n) then
+            r(3,4)=b(k+2)
+          else
+            r(3,4)=0
+          endif
+          if ((k+3).le.n) then
+           r(3,5)=d(k+3)
+          else 
+           r(3,5)=0
+          endif
+        enddo ! k
+      endsubroutine 
diff --git a/source/unres/quindibisectok.F90 b/source/unres/quindibisectok.F90
new file mode 100644 (file)
index 0000000..912d7c7
--- /dev/null
@@ -0,0 +1,68 @@
+        subroutine quindibisect2(c,b,d,n,m1,m2,eps1,relfeh,eps2,z,x) 
+        real (kind=8) :: eps1,eps2,relfeh,xmin,xmax,h,x1,xu,xo
+        real (kind=8),dimension(n):: c,b,d,x
+        real (kind=8) ,dimension(:):: wu(m1:m2)
+        integer :: i,n,m1,m2,z,a,k
+        integer maxit
+        parameter (maxit=1000)
+        
+        d(1)=0
+        d(2)=0
+        b(1)=0
+        xmin = c(n)-dabs(b(n))-dabs(d(n))
+        xmax = c(n) + dabs(b(n)) + dabs(d(n))
+        h=dabs(b(n-1))+dabs(d(n-1))+dabs(b(n))
+        if ((c(n-1)+h).gt.xmax) xmax=c(n-1)+h
+        if ((c(n-1)-h).lt.xmin) xmin=c(n-1)-h
+        do i=n-2,1,-1
+            h= dabs(b(i))+dabs(d(i))+dabs(b(i+1))+dabs(d(i+2))
+            if ((c(i)+h).gt.xmax) xmax=c(i)+h
+            if ((c(i)-h).lt.xmin) xmin=c(i)-h
+        enddo
+        if ((xmin+xmax).gt.0) then 
+           eps2=relfeh*xmax
+        else 
+           eps2=-relfeh*xmin
+        endif
+        if (eps1.lt.0) eps1=eps2
+        eps2= 0.5*eps1+7*eps2
+        xo=xmax
+        do i=m1,m2 
+          x(i)=xmax
+          wu(i)=xmin
+        enddo
+        z=0
+        do k=m2,m1,-1
+          xu=xmin
+!comment loop for the kth value
+          do i=k,m1,-1
+            if (xu.lt.wu(i)) then 
+              xu=wu(i) 
+              go to 10
+            endif
+          enddo 
+       10 if (xo.gt.x(k)) xo=x(k)
+          x1=(xu+xo)/2
+          do while  (xo-xu.gt.2*relfeh*(dabs(xu)+dabs(xo)+eps1))
+             if (z.gt.maxit) exit
+             z=z+1
+             call quindet2(c,b,d,n,x1,a)
+             a=n-a
+             if (a.lt.k) then
+                if (a.lt.m1) then 
+                   wu(m1)=x1
+                   xu=x1
+                else
+                  wu(a+1)=x1
+                  xu=x1
+                  if (x(a).gt.x1) x(a)=x1
+                endif
+             else
+                xo=x1
+             endif
+             x1=(xu+xo)/2
+          enddo
+          x(k)=(xo+xu)/2
+        enddo ! k
+       endsubroutine
+      
index 89feccc..8700369 100644 (file)
       use names
       use energy_data  !include 'COMMON.INTERACT''COMMON.CONTACTS''COMMON.VAR''COMMON.FFIELD' 'COMMON.SBRIDGE'
       use geometry_data        !include 'COMMON.GEO''COMMON.CHAIN'
- !     use REMD        !include 'COMMON.REMD'
+!     use REMD         !include 'COMMON.REMD'
 !     use MD           !include 'COMMON.MD'
       use regularize_
       use compare
 
       use io_units     !include 'COMMON.IOUNITS'
       use CSA
-
       implicit none
 #ifdef MPI
       include "mpif.h"
       real(kind=8) :: time00
       integer :: iorder,i,j,nres2,ierr,ierror
       nres2=2*nres
+#ifndef FIVEDIAG
       if(.not.allocated(Gcopy)) allocate(Gcopy(nres2,nres2))
 ! common.MD
       if(.not.allocated(Ginv)) allocate(Ginv(nres2,nres2)) !(maxres2,maxres2)
+#endif
 !      common /mdpmpi/
       if(.not.allocated(ng_start)) allocate(ng_start(0:nfgtasks-1))
       if(.not.allocated(ng_counts)) allocate(ng_counts(0:nfgtasks-1))
       if(.not.allocated(nginv_counts)) allocate(nginv_counts(0:nfgtasks-1)) !(0:MaxProcs-1)
       if(.not.allocated(nginv_start)) allocate(nginv_start(0:nfgtasks)) !(0:MaxProcs)
-
+#ifndef FIVEDIAG
       if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2)) !maxres2=2*maxres
-
+#endif
 ! Workers wait for variables and NF, and NFL from the boss
       iorder=0
       do while (iorder.ge.0)