X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2FMD_A-MTS.F;h=90c32cf2a5657a6fc19a82fe6df10e20ceeac98d;hb=83a668a2587e4624d860c7d1de07dd1634b92f54;hp=2635a2f9fc586846b39e11946fa2a0af81f69fbb;hpb=a8d848a53132a158ba39f3aa1cb20cb382d18f84;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 2635a2f..90c32cf 100644 --- a/source/unres/src_MD-M/MD_A-MTS.F +++ b/source/unres/src_MD-M/MD_A-MTS.F @@ -36,6 +36,7 @@ c------------------------------------------------ common /gucio/ cm integer itime logical ovrtim + integer nharp,iharp(4,maxres/3) c #ifdef MPI if (ilen(tmpdir).gt.0) @@ -51,6 +52,33 @@ c t_enegrad=0.0d0 t_sdsetup=0.0d0 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started" + write (iout,'(/a)') + & "Cartesian coordinates of the initial structure" + write (iout,'(a,3(3x,a5),5x,3(3x,a5))') + & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)" + do ires=1,nres + write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') + & restyp(itype(ires)),ires,(c(j,ires),j=1,3), + & (c(j,ires+nres),j=1,3) + enddo + write (iout,'(/a)') + & "Initial dC vectors of the chain" + write (iout,'(a,3(3x,a5),5x,3(3x,a5))') + & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)" + do ires=1,nres + write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') + & restyp(itype(ires)),ires,(dc(j,ires),j=1,3), + & (dc(j,ires+nres),j=1,3) + enddo + write (iout,'(/a)') + & "Initial dC_norm vectors of the chain" + write (iout,'(a,3(3x,a5),5x,3(3x,a5))') + & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)" + do ires=1,nres + write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') + & restyp(itype(ires)),ires,(dc_norm(j,ires),j=1,3), + & (dc_norm(j,ires+nres),j=1,3) + enddo #ifdef MPI tt0=MPI_Wtime() #else @@ -1443,7 +1471,7 @@ c Set up the initial conditions of a MD simulation include 'COMMON.NAMES' include 'COMMON.REMD' real*8 energia_long(0:n_ene), - & energia_short(0:n_ene),vcm(3),incr(3) + & energia_short(0:n_ene),energia(0:n_ene),vcm(3),incr(3) double precision cm(3),L(3),xv,sigv,lowb,highb double precision varia(maxvar) character*256 qstr @@ -1461,7 +1489,7 @@ c if the friction coefficients do not depend on surface area stdforcp(i)=stdfp*dsqrt(gamp) enddo do i=nnt,nct - stdforcsc(i)=stdfsc(iabs(itype(i))) + if (itype(i).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i))) & *dsqrt(gamsc(iabs(itype(i)))) enddo endif @@ -1471,9 +1499,15 @@ c Open the pdb file for snapshotshots if (ilen(tmpdir).gt.0) & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"// & liczba(:ilen(liczba))//".pdb") +#if defined(AIX) || defined(PGI) open(ipdb, & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) - & //".pdb") + & //".pdb", position='append') +#else + open(ipdb, + & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba)) + & //".pdb", access='append') +#endif else #ifdef NOXDR if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file)) @@ -1609,10 +1643,38 @@ c Removing the velocity of the center of mass write (iout,*) (vcm(j),j=1,3) endif if (.not.rest) then - call chainbuild - if(iranconf.ne.0) 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) + 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) + 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 + call etotal(energia(0)) + endif +c call chainbuild_cart + write (iout,*) "Initial energies" + call enerprint(energia(0)) + write (iout,*) "PREMINIM ",preminim + if(iranconf.ne.0 .or. preminim) then if (overlapsc) then - print *, 'Calling OVERLAP_SC' + write (iout,*) 'Calling OVERLAP_SC' call overlap_sc(fail) endif @@ -1632,8 +1694,87 @@ c Removing the velocity of the center of mass call minimize(etot,varia,iretcode,nfun) call var_to_geom(nvar,varia) endif + if(me.eq.king.or..not.out1file) then + write(iout,*) "Minimized energy is",etot + write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun + call etotal(potEcomp) + call enerprint(potEcomp) + endif + endif + if(iranconf.ne.0) then +c 8/22/17 AL Loop to produce a low-energy random conformation + DO iranmin=1,10 + if (overlapsc) then + if(me.eq.king.or..not.out1file) + & write (iout,*) 'Calling OVERLAP_SC' + call overlap_sc(fail) + endif + + 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) + 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.1.0d4) 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) +#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 endif endif call chainbuild_cart @@ -1671,11 +1812,11 @@ c Removing the velocity of the center of mass & "Time step reduced to",d_time, & " because of too large initial acceleration." endif -C if(me.eq.king.or..not.out1file)then -C write(iout,*) "Potential energy and its components" -C call enerprint(potEcomp) + 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) -C endif + endif potE=potEcomp(0)-potEcomp(20) totE=EK+potE itime=0