X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2FMD.F90;h=e3aea93843ea0e058fbfdf7ebfe67971477527d8;hb=f2b6c10194606cb27a4c4ceb30a82e6614d7e82a;hp=a03ce09fa0e1c14eb9c7bf460914581ec49896c4;hpb=d078e021535389d9b5a48b81a90e151b8978c724;p=unres4.git diff --git a/source/unres/MD.F90 b/source/unres/MD.F90 index a03ce09..e3aea93 100644 --- a/source/unres/MD.F90 +++ b/source/unres/MD.F90 @@ -813,11 +813,11 @@ !el external ilen character(len=50) :: tytul !el common /gucio/ cm - integer :: itime,i,j,nharp + integer :: i,j,nharp integer,dimension(4,nres/3) :: iharp !(4,nres/3)(4,maxres/3) ! logical :: ovrtim real(kind=8) :: tt0,scalfac - integer :: nres2 + integer :: nres2,itime nres2=2*nres print *, "ENTER MD" ! @@ -986,6 +986,7 @@ stop #endif endif + itime_mat=itime if (ntwe.ne.0) then if (mod(itime,ntwe).eq.0) then call statout(itime) @@ -1087,6 +1088,7 @@ use comm_gucio use control, only:tcpu use control_data + use minimm, only:minim_dc #ifdef MPI include 'mpif.h' integer :: ierror,ierrcode @@ -1114,17 +1116,17 @@ integer :: count,rstcount !ilen, !el external ilen character(len=50) :: tytul - integer :: maxcount_scale = 20 + integer :: maxcount_scale = 30 !el common /gucio/ cm !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres !el common /stochcalc/ stochforcvec - integer :: itime,icount_scale,itime_scal,i,j,ifac_time - logical :: scale + integer :: icount_scale,itime_scal,i,j,ifac_time,iretcode,itime + logical :: scalel real(kind=8) :: epdrift,tt0,fac_time ! if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres)) !(MAXRES6) maxres6=6*maxres - scale=.true. + scalel=.true. icount_scale=0 if (lang.eq.1) then call sddir_precalc @@ -1141,8 +1143,9 @@ #endif endif itime_scal=0 - do while (scale) + do while (scalel) icount_scale=icount_scale+1 +! write(iout,*) "icount_scale",icount_scale,scalel if (icount_scale.gt.maxcount_scale) then write (iout,*) & "ERROR: too many attempts at scaling down the time step. ",& @@ -1204,8 +1207,20 @@ call etotal(potEcomp) ! AL 4/17/17: Reduce the steps if NaNs occurred. if (potEcomp(0).gt.0.99e18 .or. isnan(potEcomp(0)).gt.0) then - d_time=d_time/2 + call enerprint(potEcomp) + d_time=d_time/10.0 + if (icount_scale.gt.15) then + write (iout,*) "Tu jest problem",potEcomp(0),d_time +! call gen_rand_conf(1,*335) +! call minim_dc(potEcomp(0),iretcode,100) + +! call zerograd +! call etotal(potEcomp) +! write(iout,*) "needed to repara,",potEcomp + endif cycle +! 335 write(iout,*) "Failed genrand" +! cycle endif ! end change if (large.and. mod(itime,ntwe).eq.0) & @@ -1230,9 +1245,13 @@ call max_accel amax=amax/(itime_scal+1)**2 call predict_edrift(epdrift) +! write(iout,*) "amax=",amax,damax,epdrift,edriftmax,amax/(itime_scal+1) + scalel=.false. +! write (iout,*) "before enter if",scalel,icount_scale if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then +! write(iout,*) "I enter if" ! Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step - scale=.true. + scalel=.true. ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax)) & /dlog(2.0d0)+1 itime_scal=itime_scal+ifac_time @@ -1348,7 +1367,6 @@ endif #endif endif - scale=.false. endif enddo ! Calculate the kinetic and the total energy and the kinetic temperature @@ -1427,7 +1445,7 @@ !el common /gucio/ cm !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres !el common /stochcalc/ stochforcvec - integer :: itime,itt,i,j,itsplit + integer :: itt,i,j,itsplit,itime logical :: scale !el common /cipiszcze/ itt @@ -2348,7 +2366,7 @@ character(len=50) :: tytul logical :: file_exist !el common /gucio/ cm - integer :: i,j,ipos,iq,iw,nft_sc,iretcode,nfun,itime,ierr,mnum + integer :: i,j,ipos,iq,iw,nft_sc,iretcode,nfun,ierr,mnum,itime real(kind=8) :: etot,tt0 logical :: fail @@ -2555,6 +2573,7 @@ if(dccart)then print *, 'Calling MINIM_DC' call minim_dc(etot,iretcode,nfun) + call int_from_cart1(.false.) else call geom_to_var(nvar,varia) print *,'Calling MINIMIZE.' @@ -2659,6 +2678,7 @@ potE=potEcomp(0)-potEcomp(20) totE=EK+potE itime=0 + itime_mat=itime if (ntwe.ne.0) call statout(itime) if(me.eq.king.or..not.out1file) & write (iout,'(/a/3(a25,1pe14.5/))') "Initial:", &