From: Cezary Czaplewski Date: Fri, 14 Apr 2017 07:36:59 +0000 (+0200) Subject: NEWGRAD added to unres Makefile and UNRES_NEWGRAD to cmake X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=commitdiff_plain;h=a01e4f4abfd7680ce30209afd6cab5d52532a55f;p=unres4.git NEWGRAD added to unres Makefile and UNRES_NEWGRAD to cmake -DUNRES_NEWGRAD=ON for cmake commandline, default OFF --- diff --git a/CMakeLists.txt b/CMakeLists.txt index 584d1d7..de3836d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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 ) #================================= diff --git a/source/unres/CMakeLists.txt b/source/unres/CMakeLists.txt index cc0b694..2a48eda 100644 --- a/source/unres/CMakeLists.txt +++ b/source/unres/CMakeLists.txt @@ -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) #========================================= diff --git a/source/unres/MD.f90 b/source/unres/MD.f90 index c509ee1..fa7131d 100644 --- a/source/unres/MD.f90 +++ b/source/unres/MD.f90 @@ -2570,6 +2570,7 @@ ! implicit real*8 (a-h,o-z) use energy_data use random, only:anorm_distr + use MD_data ! include 'DIMENSIONS' ! include 'COMMON.CONTROL' ! include 'COMMON.VAR' @@ -2588,7 +2589,15 @@ ! include 'COMMON.NAMES' ! include 'COMMON.TIME1' real(kind=8) :: xv,sigv,lowb,highb ,Ek1 - integer :: i,j,ii,k,ind +#ifdef FIVEDIAG + real(kind=8) ,allocatable, dimension(:) :: DDU1,DDU2,DL2,DL1,xsolv,DML,rs + real(kind=8) :: sumx +#ifdef DEBUG + real(kind=8) ,allocatable, dimension(:) :: rsold + real (kind=8),allocatable,dimension(:,:) :: matold +#endif +#endif + integer :: i,j,ii,k,ind,mark,imark ! Generate random velocities from Gaussian distribution of mean 0 and std of KT/m ! First generate velocities in the eigenspace of the G matrix ! write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3 @@ -2602,22 +2611,303 @@ lowb=-5*sigv highb=5*sigv d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb) -! write (iout,*) "i",i," ii",ii," geigen",geigen(i),& -! " d_t_work_new",d_t_work_new(ii) +#ifdef DEBUG + write (iout,*) "i",i," ii",ii," geigen",geigen(i),& + " d_t_work_new",d_t_work_new(ii) +#endif enddo enddo +#ifdef DEBUG ! diagnostics -! Ek1=0.0d0 -! ii=0 -! do i=1,dimen -! do k=1,3 -! ii=ii+1 -! Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2 -! enddo -! enddo -! write (iout,*) "Ek from eigenvectors",Ek1 + Ek1=0.0d0 + ii=0 + do i=1,dimen + do k=1,3 + ii=ii+1 + Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2 + enddo + enddo + write (iout,*) "Ek from eigenvectors",Ek1 + write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb) ! end diagnostics +#endif +#ifdef FIVEDIAG ! Transform velocities to UNRES coordinate space + allocate (DL1(2*nres)) + allocate (DDU1(2*nres)) + allocate (DL2(2*nres)) + allocate (DDU2(2*nres)) + allocate (xsolv(2*nres)) + allocate (DML(2*nres)) + allocate (rs(2*nres)) +#ifdef DEBUG + allocate (rsold(2*nres)) + allocate (matold(2*nres,2*nres)) + matold=0 + matold(1,1)=DMorig(1) + matold(1,2)=DU1orig(1) + matold(1,3)=DU2orig(1) + write (*,*) DMorig(1),DU1orig(1),DU2orig(1) + + do i=2,dimen-1 + do j=1,dimen + if (i.eq.j) then + matold(i,j)=DMorig(i) + matold(i,j-1)=DU1orig(i-1) + if (j.ge.3) then + matold(i,j-2)=DU2orig(i-2) + endif + + endif + + enddo + do j=1,dimen-1 + if (i.eq.j) then + matold(i,j+1)=DU1orig(i) + + end if + enddo + do j=1,dimen-2 + if(i.eq.j) then + matold(i,j+2)=DU2orig(i) + endif + enddo + enddo + matold(dimen,dimen)=DMorig(dimen) + matold(dimen,dimen-1)=DU1orig(dimen-1) + matold(dimen,dimen-2)=DU2orig(dimen-2) + write(iout,*) "old gmatrix" + call matout(dimen,dimen,2*nres,2*nres,matold) +#endif + d_t_work=0.0d0 + do i=1,dimen +! Find the ith eigenvector of the pentadiagonal inertiq matrix + do imark=dimen,1,-1 + + do j=1,imark-1 + DML(j)=DMorig(j)-geigen(i) + enddo + do j=imark+1,dimen + DML(j-1)=DMorig(j)-geigen(i) + enddo + do j=1,imark-2 + DDU1(j)=DU1orig(j) + enddo + DDU1(imark-1)=DU2orig(imark-1) + do j=imark+1,dimen-1 + DDU1(j-1)=DU1orig(j) + enddo + do j=1,imark-3 + DDU2(j)=DU2orig(j) + enddo + DDU2(imark-2)=0.0d0 + DDU2(imark-1)=0.0d0 + do j=imark,dimen-3 + DDU2(j)=DU2orig(j+1) + enddo + do j=1,dimen-3 + DL2(j+2)=DDU2(j) + enddo + do j=1,dimen-2 + DL1(j+1)=DDU1(j) + enddo +#ifdef DEBUG + write (iout,*) "DL2,DL1,DML,DDU1,DDU2" + write(iout,'(10f10.5)') (DL2(k),k=3,dimen-1) + write(iout,'(10f10.5)') (DL1(k),k=2,dimen-1) + write(iout,'(10f10.5)') (DML(k),k=1,dimen-1) + write(iout,'(10f10.5)') (DDU1(k),k=1,dimen-2) + write(iout,'(10f10.5)') (DDU2(k),k=1,dimen-3) +#endif + rs=0.0d0 + if (imark.gt.2) rs(imark-2)=-DU2orig(imark-2) + if (imark.gt.1) rs(imark-1)=-DU1orig(imark-1) + if (imark.lt.dimen) rs(imark)=-DU1orig(imark) + if (imark.lt.dimen-1) rs(imark+1)=-DU2orig(imark) +#ifdef DEBUG + rsold=rs +#endif +! write (iout,*) "Vector rs" +! do j=1,dimen-1 +! write (iout,*) j,rs(j) +! enddo + xsolv=0.0d0 + call FDIAG(dimen-1,DL2,DL1,DML,DDU1,DDU2,rs,xsolv,mark) + + if (mark.eq.1) then + +#ifdef DEBUG +! Check solution + do j=1,imark-1 + sumx=-geigen(i)*xsolv(j) + do k=1,imark-1 + sumx=sumx+matold(j,k)*xsolv(k) + enddo + do k=imark+1,dimen + sumx=sumx+matold(j,k)*xsolv(k-1) + enddo + write(iout,'(i5,3f10.5)') j,sumx,rsold(j),sumx-rsold(j) + enddo + do j=imark+1,dimen + sumx=-geigen(i)*xsolv(j-1) + do k=1,imark-1 + sumx=sumx+matold(j,k)*xsolv(k) + enddo + do k=imark+1,dimen + sumx=sumx+matold(j,k)*xsolv(k-1) + enddo + write(iout,'(i5,3f10.5)') j-1,sumx,rsold(j-1),sumx-rsold(j-1) + enddo +! end check + write (iout,*)& + "Solution of equations system",i," for eigenvalue",geigen(i) + do j=1,dimen-1 + write(iout,'(i5,f10.5)') j,xsolv(j) + enddo +#endif + do j=dimen-1,imark,-1 + xsolv(j+1)=xsolv(j) + enddo + xsolv(imark)=1.0d0 +#ifdef DEBUG + write (iout,*) "Un-normalized eigenvector",i," for eigenvalue",geigen(i) + do j=1,dimen + write(iout,'(i5,f10.5)') j,xsolv(j) + enddo +#endif +! Normalize ith eigenvector + sumx=0.0d0 + do j=1,dimen + sumx=sumx+xsolv(j)**2 + enddo + sumx=dsqrt(sumx) + do j=1,dimen + xsolv(j)=xsolv(j)/sumx + enddo +#ifdef DEBUG + write (iout,*) "Eigenvector",i," for eigenvalue",geigen(i) + do j=1,dimen + write(iout,'(i5,3f10.5)') j,xsolv(j) + enddo +#endif +! All done at this point for eigenvector i, exit loop + exit + + endif ! mark.eq.1 + + enddo ! imark + + if (mark.ne.1) then + write (iout,*) "Unable to find eigenvector",i + endif + +! write (iout,*) "i=",i + do k=1,3 +! write (iout,*) "k=",k + do j=1,dimen + ind=(j-1)*3+k +! write(iout,*) "ind",ind," ind1",3*(i-1)+k + d_t_work(ind)=d_t_work(ind) & + +xsolv(j)*d_t_work_new(3*(i-1)+k) + enddo + enddo + enddo ! i (loop over the eigenvectors) + +#ifdef DEBUG + write (iout,*) "d_t_work" + do i=1,3*dimen + write (iout,"(i5,f10.5)") i,d_t_work(i) + enddo + Ek1=0.0d0 + ii=0 + do i=nnt,nct + if (itype(i).eq.10) then + j=ii+3 + else + j=ii+6 + endif + if (i.lt.nct) then + do k=1,3 +! write (iout,*) "k",k," ii+k",ii+k," ii+j+k",ii+j+k,"EK1",Ek1 + Ek1=Ek1+0.5d0*mp*((d_t_work(ii+j+k)+d_t_work_new(ii+k))/2)**2+& + 0.5d0*0.25d0*IP*(d_t_work(ii+j+k)-d_t_work_new(ii+k))**2 + enddo + endif + if (itype(i).ne.10) ii=ii+3 + write (iout,*) "i",i," itype",itype(i)," mass",msc(itype(i)) + write (iout,*) "ii",ii + do k=1,3 + ii=ii+1 + write (iout,*) "k",k," ii",ii,"EK1",EK1 + if (iabs(itype(i)).ne.10) Ek1=Ek1+0.5d0*Isc(iabs(itype(i)))*(d_t_work(ii)-d_t_work(ii-3))**2 + Ek1=Ek1+0.5d0*msc(iabs(itype(i)))*d_t_work(ii)**2 + enddo + write (iout,*) "i",i," ii",ii + enddo + write (iout,*) "Ek from d_t_work",Ek1 + write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb) +#endif + deallocate(DDU1) + deallocate(DDU2) + deallocate(DL2) + deallocate(DL1) + deallocate(xsolv) + deallocate(DML) + deallocate(rs) +#ifdef DEBUG + deallocate(matold) + deallocate(rsold) +#endif + ind=1 + do j=nnt,nct + do k=1,3 + d_t(k,j)=d_t_work(ind) + ind=ind+1 + enddo + if (itype(j).ne.10 .and. itype(j).ne.ntyp1) then + do k=1,3 + d_t(k,j+nres)=d_t_work(ind) + ind=ind+1 + enddo + end if + enddo +#ifdef DEBUG + write (iout,*) "Random velocities in the Calpha,SC space" + do i=nnt,nct + write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3) + enddo + do i=nnt,nct + write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3) + enddo +#endif + do j=1,3 + d_t(j,0)=d_t(j,nnt) + enddo + do i=nnt,nct + if (itype(i).eq.10) then + do j=1,3 + d_t(j,i)=d_t(j,i+1)-d_t(j,i) + enddo + else + do j=1,3 + d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i) + d_t(j,i)=d_t(j,i+1)-d_t(j,i) + enddo + end if + enddo +#ifdef DEBUG + write (iout,*)"Random velocities in the virtual-bond-vector space" + do i=nnt,nct-1 + write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3) + enddo + do i=nnt,nct + write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3) + enddo + call kinetic(Ek1) + write (iout,*) "Ek from d_t_work",Ek1 + write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb) +#endif +#else do k=0,2 do i=1,dimen ind=(i-1)*3+k+1 @@ -2649,11 +2939,13 @@ enddo endif enddo +#endif ! call kinetic(EK) ! write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",& ! 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1 ! call flush(iout) return +#undef DEBUG end subroutine random_vel !----------------------------------------------------------------------------- #ifndef LANG0 @@ -5628,9 +5920,14 @@ ! common /lagrange/ allocate(d_t(3,0:nres2),d_a(3,0:nres2),d_t_old(3,0:nres2)) !(3,0:MAXRES2) allocate(d_a_work(nres6)) !(6*MAXRES) +#ifdef FIVEDIAG + allocate(DM(nres2),DU1(nres2),DU2(nres2)) + allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2)) +#else allocate(Gmat(nres2,nres2),A(nres2,nres2)) if(.not.allocated(Ginv)) allocate(Ginv(nres2,nres2)) !in control: ergastulum allocate(Gsqrp(nres2,nres2),Gsqrm(nres2,nres2),Gvec(nres2,nres2)) !(maxres2,maxres2) +#endif allocate(Geigen(nres2)) !(maxres2) if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2) ! common /inertia/ in io_conf: parmread diff --git a/source/unres/Makefile b/source/unres/Makefile index c1e9c7a..86c700c 100644 --- a/source/unres/Makefile +++ b/source/unres/Makefile @@ -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 + + + + + + + + + diff --git a/source/unres/REMD.f90 b/source/unres/REMD.f90 index edbcc8e..1deb000 100644 --- a/source/unres/REMD.f90 +++ b/source/unres/REMD.f90 @@ -41,9 +41,12 @@ ! 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 @@ -51,6 +54,123 @@ #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 @@ -70,7 +190,7 @@ 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) @@ -98,7 +218,8 @@ enddo endif enddo - + +#endif if(lmuca) then imtime=imtime+1 if(mucadyn.gt.0) call muca_update(potE) @@ -152,6 +273,7 @@ include 'mpif.h' integer :: ierror real(kind=8) :: time00 + real(kind=8) ,allocatable, dimension(:) :: DDM,DDU1,DDU2 #endif ! include 'COMMON.SETUP' ! include 'COMMON.VAR' @@ -168,10 +290,14 @@ !#endif ! include 'COMMON.IOUNITS' ! include 'COMMON.TIME1' - logical :: lprn = .false. + logical :: lprn = .true. logical :: osob +#ifndef FIVEDIAG + real(kind=8),dimension(2*nres,2*nres) :: Bmat,matmult +#endif real(kind=8) :: dtdi real(kind=8),dimension(2*nres) :: massvec,sqreig !(maxres2) maxres2=2*maxres + real(kind=8) :: relfeh,eps1,eps2 !el real(kind=8),dimension(:),allocatable :: Ghalf !el real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2)) !el real(kind=8),dimension(2*nres,2*nres) :: Gcopy !(maxres2,maxres2) @@ -180,11 +306,18 @@ integer,dimension(6*nres) :: iwork !(maxres6) maxres6=6*maxres !el common /przechowalnia/ Gcopy,Ghalf real(kind=8) :: coeff - integer :: i,j,ind,ind1,k,ii,jj,m,m1,ii1,iti,nres2,ierr + integer :: i,j,ind,ind1,k,l,ii,jj,m,m1,ii1,iti,nres2,ierr,nind,mark + real(kind=8) :: ip4 + integer :: iz + relfeh=1.0d-14 nres2=2*nres + write (iout,*) "before FIVEDIAG" +#ifndef FIVEDIAG + write (iout,*) "ALLOCATE" if(.not.allocated(Gcopy)) allocate(Gcopy(nres2,nres2)) !(maxres2,maxres2) if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !mmaxres2=(maxres2*(maxres+1)/2) +#endif ! ! Set up the matrix of the (dC,dX)-->(C,X) transformation (A), the ! inertia matrix (Gmat) and the inverse of the inertia matrix (Ginv) @@ -194,6 +327,142 @@ 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() @@ -314,6 +583,47 @@ 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 @@ -387,7 +697,8 @@ enddo enddo endif -! deallocate(Gcopy) + deallocate(Gcopy) +#endif return end subroutine setup_MD_matrices !----------------------------------------------------------------------------- diff --git a/source/unres/cinfo.f90 b/source/unres/cinfo.f90 index d38c220..22da975 100644 --- a/source/unres/cinfo.f90 +++ b/source/unres/cinfo.f90 @@ -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 diff --git a/source/unres/data/MD_data.f90 b/source/unres/data/MD_data.f90 index bee0a24..9ec0a10 100644 --- a/source/unres/data/MD_data.f90 +++ b/source/unres/data/MD_data.f90 @@ -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 diff --git a/source/unres/energy.f90 b/source/unres/energy.f90 index fdf4576..da1a2f4 100644 --- a/source/unres/energy.f90 +++ b/source/unres/energy.f90 @@ -10228,6 +10228,257 @@ 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. @@ -10466,6 +10717,7 @@ enddo return end subroutine check_ecartint +#endif !----------------------------------------------------------------------------- subroutine check_eint ! Check the gradient of energy in internal coordinates. @@ -13839,6 +14091,32 @@ (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 index 0000000..26b245f --- /dev/null +++ b/source/unres/fdiag.f90 @@ -0,0 +1,296 @@ +![ {Systems with Five--Diagonal Matrices}*) + SUBROUTINE FDIAG (N, DL2, DL1, DM, DU1, DU2, RS, X, MARK) +! +!***************************************************************** +! * +! Solving a system of linear equations * +! A * X = RS * +! with a five-diagonal, strongly nonsingular matrix A via * +! Gauss algorithm without pivoting. * +! The matrix A is given as five N-vectors DL2, DL1, DM, DU1 * +! and DU2. The linear system has the form: * +! * +! DM(1)*X(1)+DU1(1)*X(2)+DU2(1)*X(3) = RS(1) * +! DL1(2)*X(1)+DM(2)*X(2)+DU1(2)*X(3)+DU2(2)*X(4) = RS(2) * +! * +! DL2(I)*X(I-2)+DL1(I)*X(I-1)+ * +! +DM(I)*X(I)+DU1(I)*X(I+1)+DU2(I)*X(I+2) = RS(I) * +! for I = 3, ..., N - 2, and * +! * +! DL2(N-1)*X(N-3)+DL1(N-1)*X(N-2)+ * +! +DM(N-1)*X(N-1)+DU1(N-1)+X(N) = RS(N-1) * +! DL2(N)*X(N-2)+DL1(N)*X(N-1)+DM(N)*X(N) = RS(N) * +! * +! * +! * +! INPUT PARAMETERS: * +! ================= * +! N : number of equations; N > 3 * +! DL2 : N-vector DL2(1:N); second lower co-diagonal * +! DL2(3), DL2(4), ... , DL2(N) * +! DL1 : N-vector DL1(1:N); lower co-diagonal * +! DL1(2), DL1(3), ... , DL1(N) * +! DM : N-vector DM(1:N); main diagonal * +! DM(1), DM(2), ... , DM(N) * +! DU1 : N-vector DU1(1:N); upper co-diagonal * +! DU1(1), DU1(2), ... , DU1(N-1) * +! DU2 : N-vector DU2(1:N); second upper co-diagonal * +! DU2(1), DU2(2), ... , DU2(N-2) * +! RS : N-vector RS(1:N); the right hand side of the * +! linear system * +! * +! * +! OUTPUT PARAMETERS: * +! ================== * +! DL2 :) overwritten with auxiliary vectors defining the * +! DL1 :) factorization of the cyclically tridiagonal * +! DM :) matrix A * +! DU1 :) * +! DU2 :) * +! X : N-vector X(1:N); containing the solution of the * +! the system of equations * +! MARK : error parameter * +! MARK=-1 : condition N > 3 is not satisfied * +! MARK= 0 : numerically the matrix A is not strongly * +! nonsingular * +! MARK= 1 : everything is o.k. * +! * +! NOTE: if MARK = 1, the determinant of A is given by: * +! DET A = DM(1) * DM(2) * ... * DM(N) * +! * +!----------------------------------------------------------------* +! * +! subroutines required: FDIAGP, FDIAGS, MACHPD * +! * +!***************************************************************** +! * +! author : Gisela Engeln-Muellges * +! date : 05.06.1988 * +! source : FORTRAN 77 * +! * +!***************************************************************** +! + IMPLICIT DOUBLEPRECISION (A - H, O - Z) + DOUBLEPRECISION DL1 (1:N), DL2 (1:N), DM (1:N) + DOUBLEPRECISION DU1 (1:N), DU2 (1:N), RS (1:N), X (1:N) + MARK = - 1 + IF (N.LT.4) RETURN +! +! Factor the matrix A +! + CALL FDIAGP (N, DL2, DL1, DM, DU1, DU2, MARK) +! +! if MARK = 1, update and bachsubstitute +! + IF (MARK.EQ.1) THEN + CALL FDIAGS (N, DL2, DL1, DM, DU1, DU2, RS, X) + ENDIF + RETURN + END SUBROUTINE FDIAG +! +! + SUBROUTINE FDIAGP (N, DL2, DL1, DM, DU1, DU2, MARK) +! +!***************************************************************** +! * +! Factor a five-diagonal, strongly nonsingular matrix A * +! that is defined by the five N-vectors DL2, DL1, DM, DU1 * +! and DU2, into its triangular factors L * R by applying * +! Gaussian elimination specialized for five-diagonal matrices* +! (without pivoting). * +! * +! * +! INPUT PARAMETERS: * +! ================= * +! N : number of equations; N > 3 * +! DL2 : N-vector DL2(1:N); second lower co-diagonal * +! DL2(3), DL2(4), ... , DL2(N) * +! DL1 : N-vector DL1(1:N); lower co-diagonal * +! DL1(2), DL1(3), ... , DL1(N) * +! DM : N-vector DM(1:N); main diagonal * +! DM(1), DM(2), ... , DM(N) * +! DU1 : N-vector DU1(1:N); upper co-diagonal * +! DU1(1), DU1(2), ... , DU1(N-1) * +! DU2 : N-vector DU2(1:N); second upper co-diagonal * +! DU2(1), DU2(2), ... , DU2(N-2) * +! * +! * +! OUTPUT PARAMETERS: * +! ================== * +! DL2 :) overwritten with auxiliary vectors that define * +! DL1 :) the factors of the five-diagonal matrix A; * +! DM :) the three co-diagonals of the lower triangular * +! DU1 :) matrix L are stored in the vectors DL2, DL1 and * +! DU2 :) DM. The two co-diagonals of the unit upper * +! triangular matrix R are stored in the vectors DU1 * +! and DU2, its diagonal elements each have the * +! value 1. * +! MARK : error parameter * +! MARK=-1 : condition N > 3 is violated * +! MARK= 0 : numerically the matrix is not strongly * +! nonsingular * +! MARK= 1 : everything is o.k. * +! * +!----------------------------------------------------------------* +! * +! subroutines required: MACHPD * +! * +!***************************************************************** +! * +! author : Gisela Engeln-Muellges * +! date : 05.06.1988 * +! source : FORTRAN 77 * +! * +!***************************************************************** +! + IMPLICIT DOUBLEPRECISION (A - H, O - Z) + DOUBLEPRECISION DL2 (1:N), DL1 (1:N), DM (1:N), DU1 (1:N), & + DU2 (1:N) +! +! testing whether N > 3 +! + MARK = - 1 + IF (N.LT.4) RETURN +! +! calculating the machine constant +! + FMACHP = 1.0D0 + 10 FMACHP = 0.5D0 * FMACHP + IF (MACHPD (1.0D0 + FMACHP) .EQ.1) GOTO 10 + FMACHP = FMACHP * 2.0D0 +! +! determining relative error bounds +! + EPS = 4.0D0 * FMACHP +! +! initializing the undefined vector components +! + DL2 (1) = 0.0D0 + DL2 (2) = 0.0D0 + DL1 (1) = 0.0D0 + DU1 (N) = 0.0D0 + DU2 (N - 1) = 0.0D0 + DU2 (N) = 0.0D0 +! +! factoring the matrix A while checking for strong nonsingularity +! for N=1, 2 +! + ROW = DABS (DM (1) ) + DABS (DU1 (1) ) + DABS (DU2 (1) ) + IF (ROW.EQ.0.0D0) THEN + MARK = 0 + RETURN + ENDIF + D = 1.0D0 / ROW + IF (DABS (DM (1) ) * D.LE.EPS) THEN + MARK = 0 + RETURN + ENDIF + DU1 (1) = DU1 (1) / DM (1) + DU2 (1) = DU2 (1) / DM (1) + ROW = DABS (DL1 (2) ) + DABS (DM (2) ) + DABS (DU1 (2) ) + DABS ( & + DU2 (2) ) + IF (ROW.EQ.0.0D0) THEN + MARK = 0 + RETURN + ENDIF + D = 1.0D0 / ROW + DM (2) = DM (2) - DL1 (2) * DU1 (1) + IF (DABS (DM (2) ) * D.LE.EPS) THEN + MARK = 0 + RETURN + ENDIF + DU1 (2) = (DU1 (2) - DL1 (2) * DU2 (1) ) / DM (2) + DU2 (2) = DU2 (2) / DM (2) +! +! factoring A while checking for strong nonsingularity of A +! + DO 20 I = 3, N, 1 + ROW = DABS (DL2 (I) ) + DABS (DL1 (I) ) + DABS (DM (I) ) & + + DABS (DU1 (I) ) + DABS (DU2 (I) ) + IF (ROW.EQ.0.0D0) THEN + MARK = 0 + RETURN + ENDIF + D = 1.0D0 / ROW + DL1 (I) = DL1 (I) - DL2 (I) * DU1 (I - 2) + DM (I) = DM (I) - DL2 (I) * DU2 (I - 2) - DL1 (I) * DU1 (I - 1) + IF (DABS (DM (I) ) * D.LE.EPS) THEN + MARK = 0 + RETURN + ENDIF + IF (I.LT.N) THEN + DU1 (I) = (DU1 (I) - DL1 (I) * DU2 (I - 1) ) / DM (I) + ENDIF + IF (I.LT. (N - 1) ) THEN + DU2 (I) = DU2 (I) / DM (I) + ENDIF + 20 END DO + MARK = 1 + RETURN + END SUBROUTINE FDIAGP +! +! + SUBROUTINE FDIAGS (N, DL2, DL1, DM, DU1, DU2, RS, X) +! +!***************************************************************** +! * +! Solving a linear system of equations * +! A * X = RS * +! for a five-diagonal, strongly nonsingular matrix A, once * +! the factor matrices L * R have been calculated by * +! SUBROUTINE FDIAGP. Here they are used as input arrays and * +! they are stored in the five N-vectors DL2, DL1, DM, DU1 * +! and DU2. * +! * +! * +! INPUT PARAMETERS: * +! ================= * +! N : number of equations; N > 3 * +! DL2 : N-vector DL2(1:N); ) lower triangular matrix L * +! DL1 : N-vector DL1(1:N); ) including the diagonal * +! DM : N-vector DM(1:N); ) elements * +! * +! DU1 : N-vector DU1(1:N); ) unit upper triangular matrix * +! DU2 : N-vector DU2(1:N); ) R without its unit diagonal * +! elements * +! RS : N-vector RS1(1:N); right side of the linear system * +! * +! * +! OUTPUT PARAMETERS: * +! ================== * +! X : N-vector X(1:N); the solution of the linear system * +! * +!----------------------------------------------------------------* +! * +! subroutines required: none * +! * +!***************************************************************** +! * +! author : Gisela Engeln-Muellges * +! date : 05.06.1988 * +! source : FORTRAN 77 * +! * +!***************************************************************** +! + IMPLICIT DOUBLEPRECISION (A - H, O - Z) + DOUBLEPRECISION DL2 (1:N), DL1 (1:N), DM (1:N) + DOUBLEPRECISION DU1 (1:N), DU2 (1:N), RS (1:N), X (1:N) +! +! updating +! + RS (1) = RS (1) / DM (1) + RS (2) = (RS (2) - DL1 (2) * RS (1) ) / DM (2) + DO 10 I = 3, N + RS (I) = (RS (I) - DL2 (I) * RS (I - 2) - DL1 (I) * RS (I - 1) & + ) / DM (I) + 10 END DO +! +! backsubstitution +! + X (N) = RS (N) + X (N - 1) = RS (N - 1) - DU1 (N - 1) * X (N) + DO 20 I = N - 2, 1, - 1 + X (I) = RS (I) - DU1 (I) * X (I + 1) - DU2 (I) * X (I + 2) + 20 END DO + RETURN + END SUBROUTINE FDIAGS diff --git a/source/unres/fdisy.f90 b/source/unres/fdisy.f90 new file mode 100644 index 0000000..de6f33d --- /dev/null +++ b/source/unres/fdisy.f90 @@ -0,0 +1,300 @@ +![ {Systems with Five--Diagonal Symmetric Matrices} +![ {Systems with Five--Diagonal Symmetric Matrices}*) + SUBROUTINE FDISY (N, DM, DU1, DU2, RS, X, MARK) +! +!***************************************************************** +! * +! Solving a system of linear equations * +! A * X = RS * +! for a five-diagonal, symmetric and strongly nonsingular * +! matrix A. The matrix A is given by the three N-vectors DM, * +! DU1 and DU2. The system of equations has the form : * +! * +! DM(1)*X(1) + DU1(1)*X(2) + DU2(1)*X(3) = RS(1) * +! DU1(1)*X(1) + DM(2)*X(2) + DU1(2)*X(3) + DU2(2)*X(4) = RS(2) * +! * +! DU2(I-2)*X(I-2) + DU1(I-1)*X(I-1) + DM(I)*X(I) + * +! + DU1(I)*X(I+1) + DU2(I)*X(I+2) = RS(I) * +! for I = 3, ..., N - 2, and * +! * +! DU2(N-3)*X(N-2) + DU1(N-2)*X(N-1) + DM(N-1)*X(N-1) + * +! + DU1(N-1)*X(N) = RS(N-1)* +! DU2(N-2)*X(N-2) + OD(N-1)*X(N-1) + DM(N)*X(N) = RS(N) * +! * +! * +! * +! INPUT PARAMETERS: * +! ================= * +! N : number of equations, N > 3 * +! DM : N-vector DM(1:N); main diagonal of A * +! DM(1), DM(2), ... , DM(N) * +! DU1 : N-vector DU1(1:N); co-diagonal of A * +! DU1(1), DU1(2), ... , DU1(N-1) * +! DU2 : N-vector DU2(1:N); second co-diagonal of A * +! DU2(1), DU2(2), ... , DU2(N-2) * +! RS : N-vector RS(1:N); the right hand side * +! * +! * +! OUTPUT PARAMETERS: * +! ================== * +! DM :) * +! DU1 :) overwritten with intermediate quantities * +! DU2 :) * +! RS :) * +! X : N-vector X(1:N) containing the solution vector * +! MARK : error parameter * +! MARK=-2 : condition N > 3 is not satisfied * +! MARK=-1 : A is strongly nonsingular, but not positive * +! definite * +! MARK= 0 : numerically the matrix A is not strongly * +! nonsingular * +! MARK= 1 : A is positive definite * +! * +! NOTE: If MARK = +/- 1, then the determinant of A is: * +! DET A = DM(1) * DM(2) * ... * DM(N) * +! * +!----------------------------------------------------------------* +! * +! subroutines required: FDISYP, FDISYS, MACHPD * +! * +!***************************************************************** +! * +! authors : Gisela Engeln-Muellges * +! date : 01.07.1992 * +! source : FORTRAN 77 * +! * +!***************************************************************** +! + IMPLICIT DOUBLEPRECISION (A - H, O - Z) + DOUBLEPRECISION DM (1:N), DU1 (1:N), DU2 (1:N), RS (1:N), X (1:N) + MARK = - 2 + IF (N.LT.4) RETURN +! +! Factorization of the matrix A +! + CALL FDISYP (N, DM, DU1, DU2, MARK) +! +! if MARK = +/- 1 , update and backsubstitute +! + IF (MARK.EQ.1) THEN + CALL FDISYS (N, DM, DU1, DU2, RS, X) + ENDIF + RETURN + END SUBROUTINE FDISY +! +! + SUBROUTINE FDISYP (N, DM, DU1, DU2, MARK) +! +!***************************************************************** +! * +! Factor a five-diagonal, symmetric and strongly nonsingular * +! matrix A, that is given by the three N-vectors DM, DU1 and * +! DU2, into its Cholesky factors A = R(TRANSP) * D * R by * +! applying the root-free Cholesky method for five-diagonal * +! matrices. The form of the linear system is identical with * +! the one in SUBROUTINE FDISY. * +! * +! * +! INPUT PARAMETERS: * +! ================= * +! N : number of equations, N > 3 * +! DM : N-vector DM(1:N); main diagonal of A * +! DM(1), DM(2), ... , DM(N) * +! DU1 : N-vector DU1(1:N); upper co-diagonal of A * +! DU1(1), DU1(2), ... , DU1(N-1) * +! DU2 : N-vector DU2(1:N); second upper co-diagonal of A * +! DU2(1), DU2(2), ... , DU2(N-2); * +! due to symmetry the lower co-diagonals do not need to * +! be stored separately. * +! * +! * +! OUTPUT PARAMETERS: * +! ================== * +! DM :) overwritten with auxiliary vectors containing the * +! DU1 :) Cholesky factors of A. The co-diagonals of the unit * +! DU2 :) upper tridiagonal matrix R are stored in DU1 and DU2, * +! the diagonal matrix D in DM. * +! MARK : error parameter * +! MARK=-2 : condition N > 3 is not satisfied * +! MARK=-1 : A is strongly nonsingular, but not positive * +! definite * +! MARK= 0 : numerically the matrix is not strongly * +! nonsingular * +! MARK= 1 : A is positive definite * +! * +! NOTE : If MARK = +/-1, then the inertia of A, i. e., the * +! number of positive and negative eigenvalues of A, * +! is the same as the number of positive and negative * +! numbers among the components of DM. * +! * +!----------------------------------------------------------------* +! * +! subroutines required: MACHPD * +! * +!***************************************************************** +! * +! authors : Gisela Engeln-Muellges * +! date : 01.07.1988 * +! source : FORTRAN 77 * +! * +!***************************************************************** +! + IMPLICIT DOUBLEPRECISION (A - H, O - Z) + DOUBLEPRECISION DM (1:N), DU1 (1:N), DU2 (1:N) +! +! calculating the machine constant +! + FMACHP = 1.0D0 + 10 FMACHP = 0.5D0 * FMACHP + IF (MACHPD (1.0D0 + FMACHP) .EQ.1) GOTO 10 + FMACHP = FMACHP * 2.0D0 +! +! determining the relative error bound +! + EPS = 4.0D0 * FMACHP +! +! checking for N > 3 +! + MARK = - 2 + IF (N.LT.4) RETURN + DU1 (N) = 0.0D0 + DU2 (N) = 0.0D0 + DU2 (N - 1) = 0.0D0 +! +! checking for strong nonsingularity of the matrix A for N=1 +! + ROW = DABS (DM (1) ) + DABS (DU1 (1) ) + DABS (DU2 (1) ) + IF (ROW.EQ.0.0D0) THEN + MARK = 0 + RETURN + ENDIF + D = 1.0D0 / ROW + IF (DM (1) .LT.0.0D0) THEN + MARK = - 1 + RETURN + ELSEIF (DABS (DM (1) ) * D.LE.EPS) THEN + MARK = 0 + RETURN + ENDIF +! +! factoring A while checking for strong nonsingularity +! + DUMMY = DU1 (1) + DU1 (1) = DU1 (1) / DM (1) + DUMMY1 = DU2 (1) + DU2 (1) = DU2 (1) / DM (1) + ROW = DABS (DUMMY) + DABS (DM (2) ) + DABS (DU1 (2) ) + DABS (DU2 & + (2) ) + IF (ROW.EQ.0.0D0) THEN + MARK = 0 + RETURN + ENDIF + D = 1.0D0 / ROW + DM (2) = DM (2) - DUMMY * DU1 (1) + IF (DM (2) .LT.0.0D0) THEN + MARK = - 1 + RETURN + ELSEIF (DABS (DM (2) ) .LE.EPS) THEN + MARK = 0 + RETURN + ENDIF + DUMMY = DU1 (2) + DU1 (2) = (DU1 (2) - DUMMY1 * DU1 (1) ) / DM (2) + DUMMY2 = DU2 (2) + DU2 (2) = DU2 (2) / DM (2) + DO 20 I = 3, N, 1 + ROW = DABS (DUMMY1) + DABS (DUMMY) + DABS (DM (I) ) + DABS ( & + DU1 (I) ) + DABS (DU2 (I) ) + IF (ROW.EQ.0.0D0) THEN + MARK = 0 + RETURN + ENDIF + D = 1.0D0 / ROW + DM (I) = DM (I) - DM (I - 1) * DU1 (I - 1) * DU1 (I - 1) & + - DUMMY1 * DU2 (I - 2) + IF (DM (I) .LT.0.0D0) THEN + MARK = - 1 + RETURN + ELSEIF (DABS (DM (I) ) * D.LE.EPS) THEN + MARK = 0 + RETURN + ENDIF + IF (I.LT.N) THEN + DUMMY = DU1 (I) + DU1 (I) = (DU1 (I) - DUMMY2 * DU1 (I - 1) ) / DM (I) + DUMMY1 = DUMMY2 + ENDIF + IF (I.LT.N - 1) THEN + DUMMY2 = DU2 (I) + DU2 (I) = DU2 (I) / DM (I) + ENDIF + 20 END DO + MARK = 1 + RETURN + END SUBROUTINE FDISYP +! +! + SUBROUTINE FDISYS (N, DM, DU1, DU2, RS, X) +! +!***************************************************************** +! * +! Solving a linear system of equations * +! A * X = RS * +! for a five-diagonal, symmetric and strongly nonsingular * +! matrix A, once its Cholesky factors have been calculated by * +! SUBROUTINE FDISYP. Here the factors of A are used as input * +! arrays and they are stored in the three N-vectors DM, DU1 * +! and DU2. * +! * +! * +! INPUT PARAMETER: * +! ================ * +! N : number of equations, N > 3 * +! DM : N-vector DM(1:N); diagonal matrix D * +! DU1 : N-vector DM(1:N); ) co-diagonals of the upper * +! DU2 : N-vector DM(1:N); ) triangular matrix R * +! RS : N-vector DM(1:N); the right hand side * +! * +! * +! OUTPUT PARAMETER: * +! ================= * +! X : N-vector X(1:N) containing the solution vector * +! * +!----------------------------------------------------------------* +! * +! subroutines required: none * +! * +!***************************************************************** +! * +! author : Gisela Engeln-Muellges * +! date : 29.04.1988 * +! source : FORTRAN 77 * +! * +!***************************************************************** +! + IMPLICIT DOUBLEPRECISION (A - H, O - Z) + DOUBLEPRECISION DM (1:N), DU1 (1:N), DU2 (1:N), RS (1:N), X (1:N) +! +! updating +! + DUMMY1 = RS (1) + RS (1) = DUMMY1 / DM (1) + DUMMY2 = RS (2) - DU1 (1) * DUMMY1 + RS (2) = DUMMY2 / DM (2) + DO 10 I = 3, N, 1 + DUMMY1 = RS (I) - DU1 (I - 1) * DUMMY2 - DU2 (I - 2) * DUMMY1 + RS (I) = DUMMY1 / DM (I) + DUMMY3 = DUMMY2 + DUMMY2 = DUMMY1 + DUMMY1 = DUMMY3 + 10 END DO +! +! backsubstitution +! + X (N) = RS (N) + X (N - 1) = RS (N - 1) - DU1 (N - 1) * X (N) + DO 20 I = N - 2, 1, - 1 + X (I) = RS (I) - DU1 (I) * X (I + 1) - DU2 (I) * X (I + 2) + 20 END DO + RETURN + END SUBROUTINE FDISYS diff --git a/source/unres/machpd.f90 b/source/unres/machpd.f90 new file mode 100644 index 0000000..a297c99 --- /dev/null +++ b/source/unres/machpd.f90 @@ -0,0 +1,6 @@ + INTEGER FUNCTION MACHPD (X) + DOUBLEPRECISION X + MACHPD = 0 + IF (1.0D0.LT.X) MACHPD = 1 + RETURN + END FUNCTION MACHPD diff --git a/source/unres/quindet2ok.F90 b/source/unres/quindet2ok.F90 new file mode 100644 index 0000000..23e6a7b --- /dev/null +++ b/source/unres/quindet2ok.F90 @@ -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 index 0000000..912d7c7 --- /dev/null +++ b/source/unres/quindibisectok.F90 @@ -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 + diff --git a/source/unres/unres.f90 b/source/unres/unres.f90 index 89feccc..8700369 100644 --- a/source/unres/unres.f90 +++ b/source/unres/unres.f90 @@ -382,7 +382,7 @@ 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 @@ -854,7 +854,6 @@ use io_units !include 'COMMON.IOUNITS' use CSA - implicit none #ifdef MPI include "mpif.h" @@ -941,17 +940,19 @@ 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)