X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc-HCD-5D%2FMD_A-MTS.F;h=ca52aaaa65b68275c44abeeb8d5014d6f9556d80;hb=5836ecdab5a8b95f079bbf6e07374dee3fce8a26;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..ca52aaa 100644 --- a/source/unres/src-HCD-5D/MD_A-MTS.F +++ b/source/unres/src-HCD-5D/MD_A-MTS.F @@ -1111,9 +1111,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 @@ -1645,6 +1648,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) @@ -1823,7 +1832,8 @@ c Removing the velocity of the center of mass endif call flush(iout) write (iout,*) "init_MD before initial structure REST ",rest - if (.not.rest) then + 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,64 +1911,105 @@ 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 + do while (fail .and. n_model_try.lt.constr_homology) + do + i_model=iran_num(1,constr_homology) + 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 + write (iout,*) 'starting from model ',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 +#ifdef SEARCHSC + if (me.eq.king.or..not.out1file) then + write (iout,*) "Energies after removing overlaps" + call etotal(energia(0)) + call enerprint(energia(0)) + endif ! 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 + if (n_model_try.gt.constr_homology) then + write (iout,*) + & "All models have irreparable overlaps. Trying randoms starts." + iranconf=1 + goto 122 + 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 + endif ! .not. rest call chainbuild_cart call kinetic(EK) if (tbf) then call verlet_bath - endif + endif kinetic_T=2.0d0/(dimen3*Rb)*EK if(me.eq.king.or..not.out1file)then call cartprint @@ -2330,6 +2381,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