X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2FMD.f90;h=3725902b446821db4a162c33a10183125febf9ca;hb=4fa6ed0fa1ee37552df6064fe60a73f218c101ca;hp=c509ee1c7d835cf8e8a402021e24702c8e9de7d0;hpb=299e2c41124d3fa8adba7244716515a2cc160ed1;p=unres4.git diff --git a/source/unres/MD.f90 b/source/unres/MD.f90 index c509ee1..3725902 100644 --- a/source/unres/MD.f90 +++ b/source/unres/MD.f90 @@ -98,33 +98,33 @@ !el integer maxamino,maxnuc,maxbnd !el integer maxang,maxtors,maxpi !el integer maxpib,maxpit - integer :: maxatm !=2*nres !maxres2 maxres2=2*maxres - integer,parameter :: maxval=8 - integer,parameter :: maxgrp=1000 - integer,parameter :: maxtyp=3000 - integer,parameter :: maxclass=500 - integer,parameter :: maxkey=10000 - integer,parameter :: maxrot=1000 - integer,parameter :: maxopt=1000 - integer,parameter :: maxhess=1000000 - integer :: maxlight !=8*maxatm - integer,parameter :: maxvib=1000 - integer,parameter :: maxgeo=1000 - integer,parameter :: maxcell=10000 - integer,parameter :: maxring=10000 - integer,parameter :: maxfix=10000 - integer,parameter :: maxbio=10000 - integer,parameter :: maxamino=31 - integer,parameter :: maxnuc=12 - integer :: maxbnd !=2*maxatm - integer :: maxang !=3*maxatm - integer :: maxtors !=4*maxatm - integer,parameter :: maxpi=100 - integer,parameter :: maxpib=2*maxpi - integer,parameter :: maxpit=4*maxpi +! integer :: maxatm !=2*nres !maxres2 maxres2=2*maxres +! integer,parameter :: maxval=8 +! integer,parameter :: maxgrp=1000 +! integer,parameter :: maxtyp=3000 +! integer,parameter :: maxclass=500 +! integer,parameter :: maxkey=10000 +! integer,parameter :: maxrot=1000 +! integer,parameter :: maxopt=1000 +! integer,parameter :: maxhess=1000000 +! integer :: maxlight !=8*maxatm +! integer,parameter :: maxvib=1000 +! integer,parameter :: maxgeo=1000 +! integer,parameter :: maxcell=10000 +! integer,parameter :: maxring=10000 +! integer,parameter :: maxfix=10000 +! integer,parameter :: maxbio=10000 +! integer,parameter :: maxamino=31 +! integer,parameter :: maxnuc=12 +! integer :: maxbnd !=2*maxatm +! integer :: maxang !=3*maxatm +! integer :: maxtors !=4*maxatm +! integer,parameter :: maxpi=100 +! integer,parameter :: maxpib=2*maxpi +! integer,parameter :: maxpit=4*maxpi !----------------------------------------------------------------------------- ! Maximum number of seed - integer,parameter :: max_seed=1 +! integer,parameter :: max_seed=1 !----------------------------------------------------------------------------- real(kind=8),dimension(:),allocatable :: stochforcvec !(MAXRES6) maxres6=6*maxres ! common /stochcalc/ stochforcvec @@ -548,6 +548,7 @@ potE=potEcomp(0)-potEcomp(20) call cartgrad totT=totT+d_time + totTafm=totT ! Calculate the kinetic and total energy and the kinetic temperature call kinetic(EK) #ifdef MPI @@ -656,7 +657,13 @@ integer :: i,j,k,iti real(kind=8) :: KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),& mag1,mag2,v(3) - +#ifdef DEBUG + write (iout,*) "Velocities, kietic" + do i=0,nres + write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),& + (d_t(j,i+nres),j=1,3) + enddo +#endif KEt_p=0.0d0 KEt_sc=0.0d0 ! write (iout,*) "ISC",(isc(itype(i)),i=1,nres) @@ -995,10 +1002,13 @@ if (rstcount.eq.1000.or.itime.eq.n_timestep) then open(irest2,file=rest2name,status='unknown') write(irest2,*) totT,EK,potE,totE,t_bath - do i=1,2*nres + totTafm=totT +! AL 4/17/17: Now writing d_t(0,:) too + do i=0,2*nres write (irest2,'(3e15.5)') (d_t(j,i),j=1,3) enddo - do i=1,2*nres +! AL 4/17/17: Now writing d_c(0,:) too + do i=0,2*nres write (irest2,'(3e15.5)') (dc(j,i),j=1,3) enddo close(irest2) @@ -1155,6 +1165,12 @@ ! Calculate energy and forces call zerograd call etotal(potEcomp) +! AL 4/17/17: Reduce the steps if NaNs occurred. + if (potEcomp(0).gt.0.99e20 .or. isnan(potEcomp(0)).gt.0) then + d_time=d_time/2 + cycle + endif +! end change if (large.and. mod(itime,ntwe).eq.0) & call enerprint(potEcomp) #ifdef TIMING_ENE @@ -1273,6 +1289,7 @@ endif if (rattle) call rattle2 totT=totT+d_time + totTafm=totT if (d_time.ne.d_time0) then d_time=d_time0 #ifndef LANG0 @@ -1530,6 +1547,25 @@ ! Calculate energy and forces call zerograd call etotal_short(energia_short) +! AL 4/17/17: Exit itime_split loop when energy goes infinite + if (energia_short(0).gt.0.99e20 .or. isnan(energia_short(0)) ) then + if (PRINT_AMTS_MSG) & + write (iout,*) "Infinities/NaNs in energia_short",energia_short(0),"; increasing ntime_split to",ntime_split + ntime_split=ntime_split*2 + if (ntime_split.gt.maxtime_split) then +#ifdef MPI + write (iout,*) & + "Cannot rescue the run; aborting job. Retry with a smaller time step" + call flush(iout) + call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE) +#else + write (iout,*) & + "Cannot rescue the run; terminating. Retry with a smaller time step" +#endif + endif + exit + endif +! End change if (large.and. mod(itime,ntwe).eq.0) & call enerprint(energia_short) #ifdef TIMING_ENE @@ -1572,6 +1608,8 @@ if (ntime_split.lt.maxtime_split) then scale=.true. ntime_split=ntime_split*2 +! AL 4/17/17: We should exit the itime_split loop when acceleration change is too big + exit do i=0,2*nres do j=1,3 dc_old(j,i)=dc_old0(j,i) @@ -1634,6 +1672,17 @@ #endif call zerograd call etotal_long(energia_long) + if (energia_long(0).gt.0.99e20 .or. isnan(energia_long(0))) then +#ifdef MPI + write (iout,*) & + "Infinitied/NaNs in energia_long, Aborting MPI job." + call flush(iout) + call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE) +#else + write (iout,*) "Infinitied/NaNs in energia_long, terminating." + stop +#endif + endif if (large.and. mod(itime,ntwe).eq.0) & call enerprint(energia_long) #ifdef TIMING_ENE @@ -1680,6 +1729,7 @@ potE=potEcomp(0)-potEcomp(20) ! potE=energia_short(0)+energia_long(0) totT=totT+d_time + totTafm=totT ! Calculate the kinetic and the total energy and the kinetic temperature call kinetic(EK) totE=EK+potE @@ -2357,6 +2407,7 @@ endif call random_vel totT=0.0d0 + totTafm=totT endif else ! Generate initial velocities @@ -2364,6 +2415,7 @@ write(iout,*) "Initial velocities randomly generated" call random_vel totT=0.0d0 + totTafm=totT endif ! rest2name = prefix(:ilen(prefix))//'.rst' if(me.eq.king.or..not.out1file)then @@ -2389,7 +2441,7 @@ write (iout,*) "vcm right after adjustment:" write (iout,*) (vcm(j),j=1,3) endif - if (.not.rest) then + if ((.not.rest).and.(indpdb.eq.0)) then call chainbuild if(iranconf.ne.0) then if (overlapsc) then @@ -2570,6 +2622,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 +2641,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 +2663,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 +2991,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 @@ -4143,7 +4487,7 @@ !el common /przechowalnia/ GGinv,gdc,Cmat,nbond !el common /przechowalnia/ nbond integer :: max_rattle = 5 - logical :: lprn = .true., lprn1 = .true., not_done + logical :: lprn = .false., lprn1 = .false., not_done real(kind=8) :: tol_rattle = 1.0d-5 integer :: nres2 nres2=2*nres @@ -5167,11 +5511,11 @@ logical :: omit(maxarc) ! ! include 'sizes.i' - maxatm = 2*nres !maxres2 maxres2=2*maxres - maxlight = 8*maxatm - maxbnd = 2*maxatm - maxang = 3*maxatm - maxtors = 4*maxatm +! maxatm = 2*nres !maxres2 maxres2=2*maxres +! maxlight = 8*maxatm +! maxbnd = 2*maxatm +! maxang = 3*maxatm +! maxtors = 4*maxatm ! ! zero out the surface area for the sphere of interest ! @@ -5628,9 +5972,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