Merge branch 'UCGM' of mmka.chem.univ.gda.pl:unres4 into UCGM
[unres4.git] / source / unres / MD.F90
index a3451ec..3d540f5 100644 (file)
       use comm_gucio
       use control, only:tcpu
       use control_data
+      use minimm, only:minim_dc
 #ifdef MPI
       include 'mpif.h'
       integer :: ierror,ierrcode
       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 :: itime,icount_scale,itime_scal,i,j,ifac_time,iretcode
+      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
 #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. ",&
         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
+          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) &
         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
             endif
 #endif
           endif
-          scale=.false.
         endif
       enddo
 ! Calculate the kinetic and the total energy and the kinetic temperature
       character(len=16) :: form
       integer :: IERROR,ERRCODE
 #endif
+      integer :: iranmin,itrial,itmp
 !      include 'COMMON.SETUP'
 !      include 'COMMON.CONTROL'
 !      include 'COMMON.VAR'
           if(me.eq.king.or..not.out1file) &
              write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
          endif
+         if(iranconf.ne.0) then
+!c 8/22/17 AL Loop to produce a low-energy random conformation
+          DO iranmin=1,40
+          if (overlapsc) then
+           if(me.eq.king.or..not.out1file) &
+             write (iout,*) 'Calling OVERLAP_SC'
+           call overlap_sc(fail)
+          endif !endif overlap
+
+          if (searchsc) then
+           call sc_move(2,nres-1,10,1d10,nft_sc,etot)
+           print *,'SC_move',nft_sc,etot
+           if(me.eq.king.or..not.out1file) &
+           write(iout,*) 'SC_move',nft_sc,etot
+          endif
+
+          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.'
+           call minimize(etot,varia,iretcode,nfun)
+           call var_to_geom(nvar,varia)
+          endif
+          if(me.eq.king.or..not.out1file) &
+            write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
+
+          if (isnan(etot) .or. etot.gt.4.0d6) then
+            write (iout,*) "Energy too large",etot, &
+             " trying another random conformation"
+            do itrial=1,100
+              itmp=1
+              call gen_rand_conf(itmp,*30)
+              goto 40
+   30         write (iout,*) 'Failed to generate random conformation', &
+               ', itrial=',itrial
+              write (*,*) 'Processor:',me, &
+               ' Failed to generate random conformation',&
+               ' itrial=',itrial
+              call intout
+#ifdef AIX
+              call flush_(iout)
+#else
+              call flush(iout)
+#endif
+            enddo
+            write (iout,'(a,i3,a)') 'Processor:',me, &
+             ' error in generating random conformation.'
+            write (*,'(a,i3,a)') 'Processor:',me, &
+             ' error in generating random conformation.'
+            call flush(iout)
+#ifdef MPI
+!            call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
+            call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
+#else
+            stop
+#endif
+   40       continue
+          else
+            goto 44
+          endif
+          ENDDO
+
+          write (iout,'(a,i3,a)') 'Processor:',me, &
+             ' failed to generate a low-energy random conformation.'
+            write (*,'(a,i3,a,f10.3)') 'Processor:',me, &
+             ' failed to generate a low-energy random conformation.',etot
+            call flush(iout)
+            call intout
+#ifdef MPI
+!            call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
+        call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
+#else
+            stop
+#endif
+   44     continue
+        endif
       endif      
       call chainbuild_cart
       call kinetic(EK)