X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2FMD.F90;h=52d7ec3e7190cf52ab5e8802eca0386ee4847201;hb=bc23440fbe68672d430f71f22f46b11265f003db;hp=2b412fda28607e1fa5e77d71a364796bda9eb1a8;hpb=f458de69c692cc132fdb9adfa1a130fc33b35782;p=unres4.git diff --git a/source/unres/MD.F90 b/source/unres/MD.F90 index 2b412fd..52d7ec3 100644 --- a/source/unres/MD.F90 +++ b/source/unres/MD.F90 @@ -814,7 +814,7 @@ character(len=50) :: tytul !el common /gucio/ cm integer :: i,j,nharp - integer,dimension(4,nres/3) :: iharp !(4,nres/3)(4,maxres/3) + integer,dimension(4,nres) :: iharp !(4,nres/3)(4,maxres/3) ! logical :: ovrtim real(kind=8) :: tt0,scalfac integer :: nres2,itime @@ -2058,6 +2058,7 @@ do j=1,3 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time +! write(iout,*) "adt",adt,"ads",adt2 dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt d_t(j,inres)=d_t_old(j,inres)+adt @@ -2331,6 +2332,7 @@ use minimm, only:minim_dc,minimize,sc_move use io_config, only:readrst use io, only:statout + use random, only: iran_num ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' #ifdef MP @@ -2338,7 +2340,10 @@ character(len=16) :: form integer :: IERROR,ERRCODE #endif - integer :: iranmin,itrial,itmp + integer :: iranmin,itrial,itmp,n_model_try,k, & + i_model + integer, dimension(:),allocatable :: list_model_try + integer, dimension(0:nodes-1) :: i_start_models ! include 'COMMON.SETUP' ! include 'COMMON.CONTROL' ! include 'COMMON.VAR' @@ -2356,7 +2361,7 @@ ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' ! include 'COMMON.REMD' - real(kind=8),dimension(0:n_ene) :: energia_long,energia_short + real(kind=8),dimension(0:n_ene) :: energia_long,energia_short,energia real(kind=8),dimension(3) :: vcm,incr,L real(kind=8) :: xv,sigv,lowb,highb real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres) @@ -2527,8 +2532,14 @@ write (iout,*) "vcm right after adjustment:" write (iout,*) (vcm(j),j=1,3) endif + + + +! call chainbuild + if ((.not.rest).or.(forceminim)) then if (forceminim) call chainbuild_cart + 122 continue if(iranconf.ne.0 .or.indpdb.gt.0.and..not.unres_pdb .or.preminim) then if (overlapsc) then print *, 'Calling OVERLAP_SC' @@ -2633,6 +2644,136 @@ stop #endif 44 continue + else if (preminim) then + if (start_from_model) then + 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 + 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 + call int_from_cart(.true.,.false.) + call sc_loc_geom(.false.) + dc(:,0)=c(:,1) + 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 + 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) + 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 + 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 +! 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' + 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.' + call minimize(etot,varia,iretcode,nfun) + call var_to_geom(nvar,varia) +#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)& + write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun +#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 call chainbuild_cart @@ -5033,11 +5174,14 @@ !el ginvfric(2*nres,2*nres) !maxres2=2*maxres !el common /przechowalnia/ ginvfric - logical :: lprn = .false., checkmode = .false. + logical :: lprn, checkmode integer :: i,j,ind,k,nres2,nres6,mnum nres2=2*nres nres6=6*nres - + lprn=.false. + checkmode=.false. +! if (large) lprn=.true. +! if (large) checkmode=.true. if(.not.allocated(gamvec)) allocate(gamvec(nres6)) !(MAXRES6) if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres do i=0,nres2