X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2FMD.f90;h=36124aa93d3c089a732019036b66016b812a1f89;hb=19b93d89d28ac277dfba01e4e368f858db487cef;hp=fa7131db20cd082bb95ba4ac4533a1952d8c1376;hpb=a01e4f4abfd7680ce30209afd6cab5d52532a55f;p=unres4.git diff --git a/source/unres/MD.f90 b/source/unres/MD.f90 index fa7131d..36124aa 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 @@ -656,7 +656,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 +1001,12 @@ 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 +! 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 +1163,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 @@ -1530,6 +1544,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 +1605,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 +1669,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 @@ -5459,11 +5505,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 !