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})
#======================================
# 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 )
#=================================
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
)
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
#=========================================
#========================================
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)
#=========================================
! 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
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}
${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
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
+
+
+
+
+
+
+
+
+
! 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
!-----------------------------------------------------------------------------
! 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'
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
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
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
--- /dev/null
+![ {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
--- /dev/null
+![ {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
--- /dev/null
+ 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
--- /dev/null
+ 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
--- /dev/null
+ 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
+
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)