X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;ds=sidebyside;f=source%2Funres%2Fsrc-HCD-5D%2FMD_A-MTS.F;h=fcef69ea00e76ec5639360359e1c6915f5d1d131;hb=8c073a2ffd0a50442e7fdf20780a0dd94ed74799;hp=76bec98d076e58cb287c041c50c3b0f59b4afa7d;hpb=a05e8996b6fb955dc21aef7a85db21ee73d0a9d3;p=unres.git diff --git a/source/unres/src-HCD-5D/MD_A-MTS.F b/source/unres/src-HCD-5D/MD_A-MTS.F index 76bec98..fcef69e 100644 --- a/source/unres/src-HCD-5D/MD_A-MTS.F +++ b/source/unres/src-HCD-5D/MD_A-MTS.F @@ -166,14 +166,14 @@ c Entering the MD loop & .and. mod(itime,count_reset_vel).eq.0) then call random_vel write(iout,'(a,f20.2)') - & "Velocities reset to random values, time",totT + & "Velocities reset to random values, time",totT do i=0,2*nres do j=1,3 d_t_old(j,i)=d_t(j,i) enddo enddo endif - if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then + if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then call inertia_tensor call vcm_vel(vcm) do j=1,3 @@ -182,7 +182,7 @@ c Entering the MD loop call kinetic(EK) kinetic_T=2.0d0/(dimen3*Rb)*EK scalfac=dsqrt(T_bath/kinetic_T) - write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT + write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT do i=0,2*nres do j=1,3 d_t_old(j,i)=scalfac*d_t(j,i) @@ -209,6 +209,7 @@ c Variable time step algorithm. stop #endif endif + itime_mat=itime if (ntwe.ne.0) then if (mod(itime,ntwe).eq.0) then call statout(itime) @@ -263,10 +264,10 @@ C call check_ecartint 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 + do i=0,2*nres-1 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3) enddo - do i=1,2*nres + do i=0,2*nres-1 write (irest2,'(3e15.5)') (dc(j,i),j=1,3) enddo close(irest2) @@ -1111,9 +1112,12 @@ c Applying velocity Verlet algorithm - step 1 to coordinates d_t_new(j,0)=d_t_old(j,0)+adt2 d_t(j,0)=d_t_old(j,0)+adt enddo - do i=nnt,nct-1 + do i=nnt,nct-1 C SPYTAC ADAMA C do i=0,nres +#ifdef DEBUG + write (iout,*) "i",i," d_a_old",(d_a_old(j,i),j=1,3) +#endif do j=1,3 adt=d_a_old(j,i)*d_time adt2=0.5d0*adt @@ -1253,7 +1257,6 @@ c write (iout,*) "friction accelerations" call fivediaginv_mult(dimen,fric_work, d_af_work) c write (iout,*) "stochastic acceleratios" call fivediaginv_mult(dimen,stochforcvec, d_as_work) -c write (iout,*) "Leaving sddir_precalc" #else call ginv_mult(fric_work, d_af_work) call ginv_mult(stochforcvec, d_as_work) @@ -1263,6 +1266,7 @@ c write (iout,*) "Leaving sddir_precalc" write (iout,'(3f10.5)') (d_af_work(i),i=1,dimen3) write (iout,*) "d_as_work" write (iout,'(3f10.5)') (d_as_work(i),i=1,dimen3) + write (iout,*) "Leaving sddir_precalc" #endif return end @@ -1645,6 +1649,12 @@ c Set up the initial conditions of a MD simulation include 'COMMON.NAMES' include 'COMMON.REMD' include 'COMMON.TIME1' +#ifdef LBFGS + character*9 status + integer niter + common /lbfgstat/ status,niter,nfun +#endif + integer n_model_try,list_model_try(max_template),k double precision tt0 real*8 energia_long(0:n_ene), & energia_short(0:n_ene),vcm(3),incr(3) @@ -1661,6 +1671,7 @@ c Set up the initial conditions of a MD simulation integer iran_num double precision etot logical fail + integer i_start_models(0:nodes-1) write (iout,*) "init_MD INDPDB",indpdb d_time0=d_time c write(iout,*) "d_time", d_time @@ -1748,10 +1759,10 @@ c statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat' if (me.eq.king) & inquire(file=mremd_rst_name,exist=file_exist) #ifdef MPI - write (*,*) me," Before broadcast: file_exist",file_exist +c 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 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) @@ -1802,10 +1813,11 @@ c rest2name = prefix(:ilen(prefix))//'.rst' 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 + if (lang.eq.0) then c Zeroing the total angular momentum of the system write(iout,*) "Calling the zero-angular & momentum subroutine" - endif call inertia_tensor c Getting the potential energy and forces and velocities and accelerations call vcm_vel(vcm) @@ -1820,10 +1832,28 @@ c Removing the velocity of the center of mass if(me.eq.king.or..not.out1file)then write (iout,*) "vcm right after adjustment:" write (iout,*) (vcm(j),j=1,3) + write (iout,*) "Initial velocities after adjustment" + 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 + call flush(iout) endif - call flush(iout) - write (iout,*) "init_MD before initial structure REST ",rest - if (.not.rest) then + endif +c write (iout,*) "init_MD before initial structure REST ",rest + if(start_from_model .and. (me.eq.king .or. .not. out1file)) + & write(iout,*) 'START_FROM_MODELS is ON' + if(start_from_model .and. rest) then + if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then + write(iout,*) + & 'START_FROM_MODELS is OFF because the run is restarted' + write(iout,*) 'Remove restart keyword from input' + endif + endif +c write (iout,*) "rest ",rest," start_from_model",start_from_model, +c & " nmodel_start",nmodel_start," preminim",preminim + if (.not.rest) then + 122 continue if (iranconf.ne.0) then c 8/22/17 AL Loop to produce a low-energy random conformation do iranmin=1,10 @@ -1901,65 +1931,142 @@ c 8/22/17 AL Loop to produce a low-energy random conformation 44 continue else if (preminim) then if (start_from_model) then - i_model=iran_num(1,constr_homology) - write (iout,*) 'starting from model ',i_model - do i=1,2*nres - do j=1,3 - c(j,i)=chomo(j,i,i_model) + n_model_try=0 + fail=.true. + list_model_try=0 + do while (fail .and. n_model_try.lt.nmodel_start) + write (iout,*) "n_model_try",n_model_try + do + i_model=iran_num(1,nmodel_start) + do k=1,n_model_try + if (i_model.eq.list_model_try(k)) exit + enddo + if (k.gt.n_model_try) exit enddo - enddo - call int_from_cart(.true.,.false.) - call sc_loc_geom(.false.) - do i=1,nres-1 - do j=1,3 - dc(j,i)=c(j,i+1)-c(j,i) - dc_norm(j,i)=dc(j,i)*vbld_inv(i+1) + n_model_try=n_model_try+1 + list_model_try(n_model_try)=i_model + if (me.eq.king .or. .not. out1file) + & write (iout,*) 'Trying to start from model ', + & pdbfiles_chomo(i_model)(:ilen(pdbfiles_chomo(i_model))) + do i=1,2*nres + do j=1,3 + c(j,i)=chomo(j,i,i_model) + enddo enddo - enddo - do i=2,nres-1 - do j=1,3 - dc(j,i+nres)=c(j,i+nres)-c(j,i) - dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres) + call int_from_cart(.true.,.false.) + call sc_loc_geom(.false.) + do i=1,nres-1 + do j=1,3 + dc(j,i)=c(j,i+1)-c(j,i) + dc_norm(j,i)=dc(j,i)*vbld_inv(i+1) + enddo enddo - enddo - endif + do i=2,nres-1 + do j=1,3 + dc(j,i+nres)=c(j,i+nres)-c(j,i) + dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres) + enddo + enddo + if (me.eq.king.or..not.out1file) then + write (iout,*) "Energies before removing overlaps" + call etotal(energia(0)) + call enerprint(energia(0)) + endif ! Remove SC overlaps if requested - if (overlapsc) then - write (iout,*) 'Calling OVERLAP_SC' - call overlap_sc(fail) - endif + if (overlapsc) then + write (iout,*) 'Calling OVERLAP_SC' + call overlap_sc(fail) + if (fail) then + write (iout,*) + & "Failed to remove overlap from model",i_model + cycle + endif + endif + if (me.eq.king.or..not.out1file) then + write (iout,*) "Energies after removing overlaps" + call etotal(energia(0)) + call enerprint(energia(0)) + endif +#ifdef SEARCHSC ! Search for better SC rotamers if requested - 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 + 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 + call etotal(energia(0)) +#endif + enddo + call MPI_Gather(i_model,1,MPI_INTEGER,i_start_models(0), + & 1,MPI_INTEGER,king,CG_COMM,IERROR) + if (n_model_try.gt.nmodel_start .and. + & (me.eq.king .or. out1file)) then + write (iout,*) + & "All models have irreparable overlaps. Trying randoms starts." + iranconf=1 + i_model=nmodel_start+1 + goto 122 + endif + else +! Remove SC overlaps if requested + if (overlapsc) then + write (iout,*) 'Calling OVERLAP_SC' + call overlap_sc(fail) + if (fail) then + write (iout,*) + & "Failed to remove overlap" + endif + endif + if (me.eq.king.or..not.out1file) then + write (iout,*) "Energies after removing overlaps" + call etotal(energia(0)) + call enerprint(energia(0)) + endif endif - call etotal(energia(0)) C 8/22/17 AL Minimize initial structure if (dccart) then if (me.eq.king.or..not.out1file) write(iout,*) - & 'Minimizing initial PDB structure: Calling MINIM_DC' + & 'Minimizing initial PDB structure: Calling MINIM_DC' call minim_dc(etot,iretcode,nfun) else call geom_to_var(nvar,varia) if(me.eq.king.or..not.out1file) write (iout,*) - & 'Minimizing initial PDB structure: Calling MINIMIZE.' + & 'Minimizing initial PDB structure: Calling MINIMIZE.' call minimize(etot,varia,iretcode,nfun) call var_to_geom(nvar,varia) - endif - if (me.eq.king.or..not.out1file) +#ifdef LBFGS + if (me.eq.king.or..not.out1file) + & write(iout,*) 'LBFGS return code is ',status,' eval ',nfun + if(me.eq.king.or..not.out1file) + & write(iout,*) 'LBFGS return code is ',status,' eval ',nfun +#else + if (me.eq.king.or..not.out1file) & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun - if(me.eq.king.or..not.out1file) + if(me.eq.king.or..not.out1file) & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun +#endif + endif endif - endif + if (nmodel_start.gt.0 .and. me.eq.king) then + write (iout,'(a)') "Task Starting model" + do i=0,nodes-1 + if (i_start_models(i).gt.nmodel_start) then + write (iout,'(i4,2x,a)') i,"RANDOM STRUCTURE" + else + write(iout,'(i4,2x,a)')i,pdbfiles_chomo(i_start_models(i)) + & (:ilen(pdbfiles_chomo(i_start_models(i)))) + endif + enddo + endif + endif ! .not. rest call chainbuild_cart call kinetic(EK) if (tbf) then call verlet_bath - endif + endif kinetic_T=2.0d0/(dimen3*Rb)*EK + write (iout,*) "Initial kinetic energy",EK," kinetic T",kinetic_T if(me.eq.king.or..not.out1file)then call cartprint call intout @@ -1999,6 +2106,7 @@ c write(iout,*) (potEcomp(i),i=0,n_ene) c write (iout,*) "PotE-homology",potE 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:", @@ -2136,7 +2244,7 @@ c----------------------------------------------------------- double precision xv,sigv,lowb,highb,vec_afm(3),Ek1,Ek2,Ek3,aux integer i,ii,j,k,l,ind double precision anorm_distr - logical lprn /.true./ + logical lprn /.false./ #ifdef FIVEDIAG integer ichain,n,innt,inct,ibeg,ierr double precision work(8*maxres6) @@ -2300,6 +2408,9 @@ c write (iout,*) "ichain",ichain," innt",innt," inct",inct ! #define WLOS #ifdef WLOS + if (nnt.eq.1) then + d_t(:,0)=d_t(:,1) + endif do i=1,nres if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then do j=1,3 @@ -2312,11 +2423,17 @@ c write (iout,*) "ichain",ichain," innt",innt," inct",inct enddo end if enddo + d_t(:,nres)=0.0d0 d_t(:,nct)=0.0d0 + d_t(:,2*nres)=0.0d0 + if (nnt.gt.1) then + d_t(:,0)=d_t(:,1) + d_t(:,1)=0.0d0 + endif c d_a(:,0)=d_a(:,1) c d_a(:,1)=0.0d0 c write (iout,*) "Shifting accelerations" - do ichain=1,nchain + do ichain=2,nchain c write (iout,*) "ichain",chain_border1(1,ichain)-1, c & chain_border1(1,ichain) d_t(:,chain_border1(1,ichain)-1)=d_t(:,chain_border1(1,ichain)) @@ -2330,6 +2447,13 @@ c & chain_border(2,ichain-1) & d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1)) d_t(:,chain_border(2,ichain-1))=0.0d0 enddo + do ichain=2,nchain + write (iout,*) "chain",ichain,chain_border1(1,ichain)-1, + & chain_border(2,ichain-1) + d_t(:,chain_border1(1,ichain)-1)= + & d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1)) + d_t(:,chain_border(2,ichain-1))=0.0d0 + enddo #else ibeg=0 c do j=1,3 @@ -2739,8 +2863,8 @@ c enddo dc(j,0)=dc_work(j) d_t(j,0)=d_t_work(j) enddo - ind=3 - do i=nnt,nct-1 + ind=3 + do i=nnt,nct-1 do j=1,3 dc(j,i)=dc_work(ind+j) d_t(j,i)=d_t_work(ind+j)