X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2FMD_A-MTS.F;h=acd47905530f375b218bf907a5d515df80ac5f58;hb=34d3ad3987785642be58fb2f26557d3314215577;hp=6c6fb143c6ddb18cdde3ba4dac20bbcbc9411f9a;hpb=f690e8b70bab14132839afebf080d4a28363b226;p=unres.git diff --git a/source/unres/src_MD-M/MD_A-MTS.F b/source/unres/src_MD-M/MD_A-MTS.F index 6c6fb14..acd4790 100644 --- a/source/unres/src_MD-M/MD_A-MTS.F +++ b/source/unres/src_MD-M/MD_A-MTS.F @@ -196,7 +196,11 @@ c Variable time step algorithm. #endif endif if (ntwe.ne.0) then - if (mod(itime,ntwe).eq.0) call statout(itime) + if (mod(itime,ntwe).eq.0) then + call statout(itime) +C call enerprint(potEcomp) +C print *,itime,'AFM',Eafmforc,etot + endif #ifdef VOUT do j=1,3 v_work(j)=d_t(j,0) @@ -209,7 +213,7 @@ c Variable time step algorithm. enddo enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.21) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then do j=1,3 ind=ind+1 v_work(ind)=d_t(j,i+nres) @@ -230,6 +234,9 @@ c Variable time step algorithm. #endif endif if (mod(itime,ntwx).eq.0) then +c write(iout,*) 'time=',itime +C call check_ecartint + call returnbox write (tytul,'("time",f8.2)') totT if(mdpdb) then call hairpin(.true.,nharp,iharp) @@ -395,6 +402,12 @@ c Build the chain from the newly calculated coordinates c 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 @@ -513,6 +526,8 @@ c Second step of the velocity Verlet algorithm endif if (rattle) call rattle2 totT=totT+d_time + totTafm=totT +C print *,totTafm,"TU?" if (d_time.ne.d_time0) then d_time=d_time0 #ifndef LANG0 @@ -596,6 +611,8 @@ c------------------------------------------------------------------------------- include 'COMMON.IOUNITS' include 'COMMON.NAMES' include 'COMMON.TIME1' + logical lprint_short + common /shortcheck/ lprint_short double precision energia_short(0:n_ene), & energia_long(0:n_ene) double precision cm(3),L(3),vcm(3),incr(3) @@ -618,7 +635,8 @@ c------------------------------------------------------------------------------- write (iout,*) "***************** RESPA itime",itime write (iout,*) "Cartesian and internal coordinates: step 0" c call cartprint - call pdbout(0.0d0,"cipiszcze",iout) + call pdbout(0.0d0, + & "cipiszcze ",iout) call intout write (iout,*) "Accelerations from long-range forces" do i=0,nres @@ -742,7 +760,8 @@ c Build the chain from the newly calculated coordinates if (large.and. mod(itime,ntwe).eq.0) then write (iout,*) "***** ITSPLIT",itsplit write (iout,*) "Cartesian and internal coordinates: step 1" - call pdbout(0.0d0,"cipiszcze",iout) + call pdbout(0.0d0, + & "cipiszcze ",iout) c call cartprint call intout write (iout,*) "Velocities, step 1" @@ -758,8 +777,30 @@ c call cartprint tt0 = tcpu() #endif c Calculate energy and forces +c if (large.and. mod(itime,ntwe).eq.0) lprint_short=.true. 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 @@ -864,6 +905,22 @@ c Compute long-range forces #endif call zerograd call etotal_long(energia_long) +! AL 4/17/2017 Handling NaNs + 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 enerprint(energia_long) + call flush(iout) + call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE) +#else + write (iout,*) "Infinitied/NaNs in energia_long, terminating." + call enerprint(energia_long) + stop +#endif + endif +! end change +c lprint_short=.false. if (large.and. mod(itime,ntwe).eq.0) & call enerprint(energia_long) #ifdef TIMING_ENE @@ -886,7 +943,8 @@ c Compute accelerations from long-range forces write (iout,*) "energia_long",energia_long(0) write (iout,*) "Cartesian and internal coordinates: step 2" c call cartprint - call pdbout(0.0d0,"cipiszcze",iout) + call pdbout(0.0d0, + & cipiszcze ,iout) call intout write (iout,*) "Accelerations from long-range forces" do i=0,nres @@ -908,8 +966,15 @@ c Compute the complete potential energy potEcomp(i)=energia_short(i)+energia_long(i) enddo potE=potEcomp(0)-potEcomp(20) + if (ntwe.ne.0) then + if (large.and. mod(itime,ntwe).eq.0) then + call enerprint(potEcomp) + write (iout,*) "potE",potD + endif + endif c potE=energia_short(0)+energia_long(0) totT=totT+d_time + totTafm=totT c Calculate the kinetic and the total energy and the kinetic temperature call kinetic(EK) totE=EK+potE @@ -955,7 +1020,7 @@ c forces). enddo enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.21) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then inres=i+nres do j=1,3 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time @@ -996,6 +1061,8 @@ c Applying velocity Verlet algorithm - step 1 to coordinates d_t(j,0)=d_t_old(j,0)+adt enddo do i=nnt,nct-1 +C SPYTAC ADAMA +C do i=0,nres do j=1,3 adt=d_a_old(j,i)*d_time adt2=0.5d0*adt @@ -1005,7 +1072,8 @@ c Applying velocity Verlet algorithm - step 1 to coordinates enddo enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.21) then +C do i=0,nres + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then inres=i+nres do j=1,3 adt=d_a_old(j,inres)*d_time @@ -1049,7 +1117,7 @@ c Step 2 of the velocity Verlet algorithm: update velocities enddo enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.21) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then inres=i+nres do j=1,3 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time @@ -1087,12 +1155,25 @@ c Applying velocity Verlet algorithm - step 1 to coordinates c c Compute friction and stochastic forces c +#ifdef MPI time00=MPI_Wtime() +#else + time00=tcpu() +#endif call friction_force +#ifdef MPI time_fric=time_fric+MPI_Wtime()-time00 time00=MPI_Wtime() +#else + time_fric=time_fric+tcpu()-time00 + time00=tcpu() +#endif call stochastic_force(stochforcvec) +#ifdef MPI time_stoch=time_stoch+MPI_Wtime()-time00 +#else + time_stoch=time_stoch+tcpu()-time00 +#endif c c Compute the acceleration due to friction forces (d_af_work) and stochastic c forces (d_as_work) @@ -1148,7 +1229,7 @@ c ind=ind+3 enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.21) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then inres=i+nres do j=1,3 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time @@ -1213,7 +1294,7 @@ c ind=ind+3 enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.21) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then inres=i+nres do j=1,3 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres) @@ -1260,6 +1341,7 @@ c accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i)) c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j)) if (dabs(accel(j)).gt.dabs(accel_old(j))) then dacc=dabs(accel(j)-accel_old(j)) +c write (iout,*) i,dacc if (dacc.gt.amax) amax=dacc endif enddo @@ -1278,7 +1360,7 @@ c accel(j)=aux(j) enddo endif do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.21) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then do j=1,3 c accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres) accel_old(j)=accel_old(j)+d_a_old(j,i+nres) @@ -1289,6 +1371,7 @@ c accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres) c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j)) if (dabs(accel(j)).gt.dabs(accel_old(j))) then dacc=dabs(accel(j)-accel_old(j)) +c write (iout,*) "side-chain",i,dacc if (dacc.gt.amax) amax=dacc endif enddo @@ -1331,7 +1414,7 @@ c write (iout,*) "back",i,j,epdriftij enddo endif c Side chains - if (itype(i).ne.10 .and. itype(i).ne.21) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then do j=1,3 epdriftij= & dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i)) @@ -1378,7 +1461,7 @@ c write(iout,*) "fact", fact enddo enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.21) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then inres=i+nres do j=1,3 d_t(j,inres)=fact*d_t(j,inres) @@ -1433,7 +1516,8 @@ c if the friction coefficients do not depend on surface area stdforcp(i)=stdfp*dsqrt(gamp) enddo do i=nnt,nct - stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(itype(i))) + stdforcsc(i)=stdfsc(iabs(itype(i))) + & *dsqrt(gamsc(iabs(itype(i)))) enddo endif c Open the pdb file for snapshotshots @@ -1507,11 +1591,13 @@ c statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat' if (restart1file) then if (me.eq.king) & inquire(file=mremd_rst_name,exist=file_exist) +#ifdef MPI write (*,*) me," Before broadcast: file_exist",file_exist call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM, & IERR) write (*,*) me," After broadcast: file_exist",file_exist c inquire(file=mremd_rst_name,exist=file_exist) +#endif if(me.eq.king.or..not.out1file) & write(iout,*) "Initial state read by master and distributed" else @@ -1529,6 +1615,7 @@ c inquire(file=mremd_rst_name,exist=file_exist) endif call rescale_weights(t_bath) else + rest=.false. if(me.eq.king.or..not.out1file)then if (restart1file) then write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)), @@ -1541,6 +1628,7 @@ c inquire(file=mremd_rst_name,exist=file_exist) endif call random_vel totT=0.0d0 + totTafm=totT endif else c Generate initial velocities @@ -1548,6 +1636,8 @@ c Generate initial velocities & write(iout,*) "Initial velocities randomly generated" call random_vel totT=0.0d0 +CtotTafm is the variable for AFM time which eclipsed during + totTafm=totT endif c rest2name = prefix(:ilen(prefix))//'.rst' if(me.eq.king.or..not.out1file)then @@ -1559,24 +1649,40 @@ c rest2name = prefix(:ilen(prefix))//'.rst' c Zeroing the total angular momentum of the system write(iout,*) "Calling the zero-angular & momentum subroutine" + call flush(iout) endif call inertia_tensor +c write (iout,*) "After inertia" +c call flush(iout) c Getting the potential energy and forces and velocities and accelerations + if(me.eq.king.or..not.out1file)then + write(iout,*) "Calling the vcm_vel" + call flush(iout) + endif call vcm_vel(vcm) -c write (iout,*) "velocity of the center of the mass:" -c write (iout,*) (vcm(j),j=1,3) + write (iout,*) "velocity of the center of the mass:" + write (iout,*) (vcm(j),j=1,3) + call flush(iout) do j=1,3 d_t(j,0)=d_t(j,0)-vcm(j) enddo c Removing the velocity of the center of mass + if(me.eq.king.or..not.out1file)then + write(iout,*) "Calling the vcm_vel" + call flush(iout) + endif call vcm_vel(vcm) if(me.eq.king.or..not.out1file)then write (iout,*) "vcm right after adjustment:" write (iout,*) (vcm(j),j=1,3) + call flush(iout) endif if (.not.rest) then - call chainbuild if(iranconf.ne.0) then +c 8/22/17 AL Loop to produce a low-energy random conformation + do iranmin=1,10 + write (iout,*) "iranmin",iranmin + call chainbuild if (overlapsc) then print *, 'Calling OVERLAP_SC' call overlap_sc(fail) @@ -1600,12 +1706,72 @@ c Removing the velocity of the center of mass endif if(me.eq.king.or..not.out1file) & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun + if (isnan(etot) .or. etot.gt.1.0d4) then + write (iout,*) "Energy too large",etot, + & " trying another random conformation" + do itrial=1,100 + itmp=1 + nrestmp=nres + call gen_rand_conf(itmp,nrestmp,*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) +#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)') 'Processor:',me, + & ' failed to generate a low-energy random conformation.' + call flush(iout) +#ifdef MPI + call MPI_Abort(mpi_comm_world,error_msg,ierrcode) +#else + stop +#endif + 44 continue + else if (indpdb.gt.0) then +C 8/22/17 AL Minimize initial PDB structure + if(dccart)then + print *, 'Calling MINIM_DC' + call minim_dc(etot,iretcode,nfun) + 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 endif endif call chainbuild_cart call kinetic(EK) if (tbf) then - call verlet_bath(EK) + call verlet_bath endif kinetic_T=2.0d0/(dimen3*Rb)*EK if(me.eq.king.or..not.out1file)then @@ -1640,7 +1806,7 @@ c Removing the velocity of the center of mass if(me.eq.king.or..not.out1file)then write(iout,*) "Potential energy and its components" call enerprint(potEcomp) -c write(iout,*) (potEcomp(i),i=0,n_ene) + write(iout,*) (potEcomp(i),i=0,n_ene) endif potE=potEcomp(0)-potEcomp(20) totE=EK+potE @@ -1770,7 +1936,7 @@ c----------------------------------------------------------- include 'COMMON.IOUNITS' include 'COMMON.NAMES' include 'COMMON.TIME1' - double precision xv,sigv,lowb,highb + double precision xv,sigv,lowb,highb,vec_afm(3) c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m c First generate velocities in the eigenspace of the G matrix c write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3 @@ -1784,10 +1950,27 @@ c call flush(iout) lowb=-5*sigv highb=5*sigv d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb) + c write (iout,*) "i",i," ii",ii," geigen",geigen(i), c & " d_t_work_new",d_t_work_new(ii) enddo enddo +C if (SELFGUIDE.gt.0) then +C distance=0.0 +C do j=1,3 +C vec_afm(j)=c(j,afmend)-c(j,afmbeg) +C distance=distance+vec_afm(j)**2 +C enddo +C distance=dsqrt(distance) +C do j=1,3 +C d_t_work_new(j+(afmbeg-1)*3)=-velAFMconst*vec_afm(j)/distance +C d_t_work_new(j+(afmend-1)*3)=velAFMconst*vec_afm(j)/distance +C write(iout,*) "myvel",d_t_work_new(j+(afmbeg-1)*3), +C & d_t_work_new(j+(afmend-1)*3) +C enddo + +C endif + c diagnostics c Ek1=0.0d0 c ii=0 @@ -1824,7 +2007,7 @@ c Transfer to the d_t vector enddo enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.21) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then do j=1,3 ind=ind+1 d_t(j,i+nres)=d_t_work(ind) @@ -2053,7 +2236,7 @@ c enddo ind=ind+3 enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.21) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then do j=1,3 dc_work(ind+j)=dc_old(j,i+nres) d_t_work(ind+j)=d_t_old(j,i+nres) @@ -2162,7 +2345,7 @@ c ind=ind+3 enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.21) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then inres=i+nres do j=1,3 d_t(j,inres)=d_t_work(ind+j) @@ -2319,7 +2502,7 @@ c enddo ind=ind+3 enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.21) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then do j=1,3 dc_work(ind+j)=dc_old(j,i+nres) d_t_work(ind+j)=d_t_old(j,i+nres) @@ -2368,7 +2551,7 @@ c enddo ind=ind+3 enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.21) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then inres=i+nres do j=1,3 dc(j,inres)=dc_work(ind+j) @@ -2429,7 +2612,7 @@ c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j) ind=ind+3 enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.21) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then inres=i+nres do j=1,3 d_t(j,inres)=d_t_work(ind+j)