X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc-HCD-5D%2FMD_A-MTS.F;h=d82cf17827596cd63f8b82bd27f66ed2c6c77e57;hb=a9f1b123976804db7a99d09f2f43296d9fd1cf0d;hp=a8efa206e66695419a1609373efd3568fbe724c5;hpb=020e579626d686ec20ecd9f0cc4c8313f474e152;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 a8efa20..d82cf17 100644 --- a/source/unres/src-HCD-5D/MD_A-MTS.F +++ b/source/unres/src-HCD-5D/MD_A-MTS.F @@ -2,7 +2,7 @@ c------------------------------------------------ c The driver for molecular dynamics subroutines c------------------------------------------------ - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' #ifdef MPI include "mpif.h" @@ -12,11 +12,20 @@ c------------------------------------------------ include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' +#ifdef FIVEDIAG + include 'COMMON.LAGRANGE.5diag' +#else + include 'COMMON.LAGRANGE' +#endif #ifndef LANG0 include 'COMMON.LANGEVIN' #else +#ifdef FIVEDIAG + include 'COMMON.LANGEVIN.lang0.5diag' +#else include 'COMMON.LANGEVIN.lang0' #endif +#endif include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' @@ -36,6 +45,10 @@ c------------------------------------------------ common /gucio/ cm integer itime logical ovrtim + integer i,j,icount_scale,itime_scal + integer nharp,iharp(4,maxres/3) + double precision scalfac + double precision tt0 c #ifdef MPI if (ilen(tmpdir).gt.0) @@ -45,6 +58,7 @@ c if (ilen(tmpdir).gt.0) & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst') #endif + write (iout,*) "MD lang",lang t_MDsetup=0.0d0 t_langsetup=0.0d0 t_MD=0.0d0 @@ -152,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 @@ -168,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) @@ -195,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) @@ -288,21 +303,30 @@ c------------------------------------------------------------------------------- c Perform a single velocity Verlet step; the time step can be rescaled if c increments in accelerations exceed the threshold c------------------------------------------------------------------------------- - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' #ifdef MPI include 'mpif.h' - integer ierror,ierrcode + integer ierror,ierrcode,errcode #endif include 'COMMON.SETUP' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' +#ifdef FIVEDIAG + include 'COMMON.LAGRANGE.5diag' +#else + include 'COMMON.LAGRANGE' +#endif #ifndef LANG0 include 'COMMON.LANGEVIN' #else +#ifdef FIVEDIAG + include 'COMMON.LANGEVIN.lang0.5diag' +#else include 'COMMON.LANGEVIN.lang0' #endif +#endif include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' @@ -321,8 +345,10 @@ c------------------------------------------------------------------------------- common /gucio/ cm double precision stochforcvec(MAXRES6) common /stochcalc/ stochforcvec - integer itime + integer itime,icount_scale,itime_scal,ifac_time,i,j,itt logical scale + double precision epdrift,fac_time + double precision tt0 c scale=.true. icount_scale=0 @@ -403,7 +429,7 @@ 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))) then + if (potEcomp(0).gt.0.99e20 .or. isnan(potEcomp(0)).gt.0) then d_time=d_time/2 cycle endif @@ -588,7 +614,7 @@ c------------------------------------------------------------------------------- c------------------------------------------------------------------------------- c Perform a single RESPA step. c------------------------------------------------------------------------------- - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' #ifdef MPI include 'mpif.h' @@ -598,11 +624,20 @@ c------------------------------------------------------------------------------- include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' +#ifdef FIVEDIAG + include 'COMMON.LAGRANGE.5diag' +#else + include 'COMMON.LAGRANGE' +#endif #ifndef LANG0 include 'COMMON.LANGEVIN' #else +#ifdef FIVEDIAG + include 'COMMON.LANGEVIN.lang0.5diag' +#else include 'COMMON.LANGEVIN.lang0' #endif +#endif include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' @@ -618,6 +653,8 @@ c------------------------------------------------------------------------------- double precision cm(3),L(3),vcm(3),incr(3) double precision dc_old0(3,0:maxres2),d_t_old0(3,0:maxres2), & d_a_old0(3,0:maxres2) + integer i,j + double precision fac_time logical PRINT_AMTS_MSG /.false./ integer ilen,count,rstcount external ilen @@ -628,7 +665,11 @@ c------------------------------------------------------------------------------- common /stochcalc/ stochforcvec integer itime logical scale + integer itt common /cipiszcze/ itt + integer itsplit + double precision epdrift,epdriftmax + double precision tt0 itt=itime if (ntwe.ne.0) then if (large.and. mod(itime,ntwe).eq.0) then @@ -944,7 +985,7 @@ c Compute accelerations from long-range forces write (iout,*) "Cartesian and internal coordinates: step 2" c call cartprint call pdbout(0.0d0, - & cipiszcze ,iout) + & 'cipiszcze ',iout) call intout write (iout,*) "Accelerations from long-range forces" do i=0,nres @@ -969,7 +1010,7 @@ c Compute the complete potential energy if (ntwe.ne.0) then if (large.and. mod(itime,ntwe).eq.0) then call enerprint(potEcomp) - write (iout,*) "potE",potD + write (iout,*) "potE",potE endif endif c potE=energia_short(0)+energia_long(0) @@ -999,11 +1040,16 @@ c--------------------------------------------------------------------- subroutine RESPA_vel c First and last RESPA step (incrementing velocities using long-range c forces). - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' +#ifdef FIVEDIAG + include 'COMMON.LAGRANGE.5diag' +#else + include 'COMMON.LAGRANGE' +#endif include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' @@ -1011,6 +1057,7 @@ c forces). include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' + integer i,j,inres do j=1,3 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time enddo @@ -1032,11 +1079,16 @@ c forces). c----------------------------------------------------------------- subroutine verlet1 c Applying velocity Verlet algorithm - step 1 to coordinates - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' +#ifdef FIVEDIAG + include 'COMMON.LAGRANGE.5diag' +#else + include 'COMMON.LAGRANGE' +#endif include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' @@ -1045,7 +1097,7 @@ c Applying velocity Verlet algorithm - step 1 to coordinates include 'COMMON.IOUNITS' include 'COMMON.NAMES' double precision adt,adt2 - + integer i,j,inres #ifdef DEBUG write (iout,*) "VELVERLET1 START: DC" do i=0,nres @@ -1060,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 @@ -1096,11 +1151,16 @@ C do i=0,nres c--------------------------------------------------------------------- subroutine verlet2 c Step 2 of the velocity Verlet algorithm: update velocities - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' +#ifdef FIVEDIAG + include 'COMMON.LAGRANGE.5diag' +#else + include 'COMMON.LAGRANGE' +#endif include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' @@ -1108,6 +1168,7 @@ c Step 2 of the velocity Verlet algorithm: update velocities include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' + integer i,j,inres do j=1,3 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time enddo @@ -1129,7 +1190,7 @@ c Step 2 of the velocity Verlet algorithm: update velocities c----------------------------------------------------------------- subroutine sddir_precalc c Applying velocity Verlet algorithm - step 1 to coordinates - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' #ifdef MPI include 'mpif.h' @@ -1137,11 +1198,20 @@ c Applying velocity Verlet algorithm - step 1 to coordinates include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' +#ifdef FIVEDIAG + include 'COMMON.LAGRANGE.5diag' +#else + include 'COMMON.LAGRANGE' +#endif #ifndef LANG0 include 'COMMON.LANGEVIN' #else +#ifdef FIVEDIAG + include 'COMMON.LANGEVIN.lang0.5diag' +#else include 'COMMON.LANGEVIN.lang0' #endif +#endif include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' @@ -1150,8 +1220,10 @@ c Applying velocity Verlet algorithm - step 1 to coordinates include 'COMMON.IOUNITS' include 'COMMON.NAMES' include 'COMMON.TIME1' + double precision time00 double precision stochforcvec(MAXRES6) common /stochcalc/ stochforcvec + integer i c c Compute friction and stochastic forces c @@ -1161,6 +1233,7 @@ c time00=tcpu() #endif call friction_force +c write (iout,*) "After friction_force" #ifdef MPI time_fric=time_fric+MPI_Wtime()-time00 time00=MPI_Wtime() @@ -1169,6 +1242,7 @@ c time00=tcpu() #endif call stochastic_force(stochforcvec) +c write (iout,*) "After stochastic_force" #ifdef MPI time_stoch=time_stoch+MPI_Wtime()-time00 #else @@ -1178,23 +1252,46 @@ c c Compute the acceleration due to friction forces (d_af_work) and stochastic c forces (d_as_work) c +#ifdef FIVEDIAG +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) +#else call ginv_mult(fric_work, d_af_work) call ginv_mult(stochforcvec, d_as_work) +#endif +#ifdef DEBUG + write (iout,*) "d_af_work" + 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 c--------------------------------------------------------------------- subroutine sddir_verlet1 c Applying velocity Verlet algorithm - step 1 to velocities - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' +#ifdef FIVEDIAG + include 'COMMON.LAGRANGE.5diag' +#else + include 'COMMON.LAGRANGE' +#endif #ifndef LANG0 include 'COMMON.LANGEVIN' #else +#ifdef FIVEDIAG + include 'COMMON.LANGEVIN.lang0.5diag' +#else include 'COMMON.LANGEVIN.lang0' #endif +#endif include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' @@ -1206,6 +1303,7 @@ c Revised 3/31/05 AL: correlation between random contributions to c position and velocity increments included. double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3) double precision adt,adt2 + integer i,j,ind,inres c c Add the contribution from BOTH friction and stochastic force to the c coordinates, but ONLY the contribution from the friction forces to velocities @@ -1218,7 +1316,7 @@ c d_t(j,0)=d_t_old(j,0)+adt enddo ind=3 - do i=nnt,nct-1 + do i=nnt,nct-1 do j=1,3 adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time @@ -1246,16 +1344,25 @@ c c--------------------------------------------------------------------- subroutine sddir_verlet2 c Calculating the adjusted velocities for accelerations - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' +#ifdef FIVEDIAG + include 'COMMON.LAGRANGE.5diag' +#else + include 'COMMON.LAGRANGE' +#endif #ifndef LANG0 include 'COMMON.LANGEVIN' #else +#ifdef FIVEDIAG + include 'COMMON.LANGEVIN.lang0.5diag' +#else include 'COMMON.LANGEVIN.lang0' #endif +#endif include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' @@ -1265,6 +1372,7 @@ c Calculating the adjusted velocities for accelerations include 'COMMON.NAMES' double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6) double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/ + integer i,j,inres,ind c Revised 3/31/05 AL: correlation between random contributions to c position and velocity increments included. c The correlation coefficients are calculated at low-friction limit. @@ -1276,8 +1384,11 @@ c c Compute the acceleration due to friction forces (d_af_work) and stochastic c forces (d_as_work) c +#ifdef FIVEDIAG + call fivediaginv_mult(maxres6,stochforcvec, d_as_work1) +#else call ginv_mult(stochforcvec, d_as_work1) - +#endif c c Update velocities c @@ -1312,17 +1423,23 @@ c c Find the maximum difference in the accelerations of the the sites c at the beginning and the end of the time step. c - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' +#ifdef FIVEDIAG + include 'COMMON.LAGRANGE.5diag' +#else + include 'COMMON.LAGRANGE' +#endif include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' + integer i,j double precision aux(3),accel(3),accel_old(3),dacc do j=1,3 c aux(j)=d_a(j,0)-d_a_old(j,0) @@ -1388,11 +1505,16 @@ c--------------------------------------------------------------------- c c Predict the drift of the potential energy c - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' +#ifdef FIVEDIAG + include 'COMMON.LAGRANGE.5diag' +#else + include 'COMMON.LAGRANGE' +#endif include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' @@ -1401,6 +1523,7 @@ c include 'COMMON.IOUNITS' include 'COMMON.MUCA' double precision epdrift,epdriftij + integer i,j c Drift of the potential energy epdrift=0.0d0 do i=nnt,nct @@ -1433,11 +1556,25 @@ c----------------------------------------------------------------------- c c Coupling to the thermostat by using the Berendsen algorithm c - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' +#ifdef FIVEDIAG + include 'COMMON.LAGRANGE.5diag' +#else + include 'COMMON.LAGRANGE' +#endif +#ifdef LANG0 +#ifdef FIVEDIAG + include 'COMMON.LANGEVIN.lang0.5diag' +#else + include 'COMMON.LANGEVIN.lang0' +#endif +#else + include 'COMMON.LANGEVIN' +#endif include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' @@ -1446,7 +1583,7 @@ c include 'COMMON.IOUNITS' include 'COMMON.NAMES' double precision T_half,fact -c + integer i,j ,inres T_half=2.0d0/(dimen3*Rb)*EK fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0)) c write(iout,*) "T_half", T_half @@ -1473,22 +1610,36 @@ c write(iout,*) "fact", fact c--------------------------------------------------------- subroutine init_MD c Set up the initial conditions of a MD simulation - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' #ifdef MP include 'mpif.h' character*16 form - integer IERROR,ERRCODE + integer IERROR,ERRCODE,error_msg,ierr,ierrcode #endif include 'COMMON.SETUP' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' +#ifdef FIVEDIAG + include 'COMMON.LAGRANGE.5diag' +#else + include 'COMMON.LAGRANGE' +#endif + include 'COMMON.QRESTR' #ifndef LANG0 include 'COMMON.LANGEVIN' #else +#ifdef FIVEDIAG + include 'COMMON.LANGEVIN.lang0.5diag' +#else +#ifdef FIVEDIAG + include 'COMMON.LANGEVIN.lang0.5diag' +#else include 'COMMON.LANGEVIN.lang0' #endif +#endif +#endif include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' @@ -1497,6 +1648,14 @@ c Set up the initial conditions of a MD simulation include 'COMMON.IOUNITS' 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) double precision cm(3),L(3),xv,sigv,lowb,highb @@ -1507,6 +1666,11 @@ c Set up the initial conditions of a MD simulation character*50 tytul logical file_exist common /gucio/ cm + integer i,ipos,iq,iw,j,iranmin,nft_sc,iretcode,nfun,itrial,itmp, + & i_model,itime + integer iran_num + double precision etot + logical fail write (iout,*) "init_MD INDPDB",indpdb d_time0=d_time c write(iout,*) "d_time", d_time @@ -1648,10 +1812,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) @@ -1666,10 +1831,17 @@ 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 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 @@ -1747,64 +1919,120 @@ 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 + 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 + if (n_model_try.gt.constr_homology) then + write (iout,*) + & "All models have irreparable overlaps. Trying randoms starts." + iranconf=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 + 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 @@ -1845,6 +2073,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:", @@ -1952,16 +2181,25 @@ C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array end c----------------------------------------------------------- subroutine random_vel - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' +#ifdef FIVEDIAG + include 'COMMON.LAGRANGE.5diag' +#else + include 'COMMON.LAGRANGE' +#endif #ifndef LANG0 include 'COMMON.LANGEVIN' #else +#ifdef FIVEDIAG + include 'COMMON.LANGEVIN.lang0.5diag' +#else include 'COMMON.LANGEVIN.lang0' #endif +#endif include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' @@ -1970,11 +2208,281 @@ c----------------------------------------------------------- include 'COMMON.IOUNITS' include 'COMMON.NAMES' include 'COMMON.TIME1' - double precision xv,sigv,lowb,highb,vec_afm(3) + 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 /.false./ +#ifdef FIVEDIAG + integer ichain,n,innt,inct,ibeg,ierr + double precision work(8*maxres6) + integer iwork(maxres6) + double precision Ghalf(mmaxres2_chain),Geigen(maxres2_chain), + & Gvec(maxres2_chain,maxres2_chain) + common /przechowalnia/Ghalf,Geigen,Gvec +#ifdef DEBUG + double precision inertia(maxres2_chain,maxres2_chain) +#endif 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 c call flush(iout) +#ifdef DEBUG + write (iout,*) "Random_vel, fivediag" +#endif + d_t=0.0d0 + Ek2=0.0d0 + EK=0.0d0 + Ek3=0.0d0 + do ichain=1,nchain + ind=0 + ghalf=0.0d0 + n=dimen_chain(ichain) + innt=iposd_chain(ichain) + inct=innt+n-1 +#ifdef DEBUG + write (iout,*) "Chain",ichain," n",n," start",innt + do i=innt,inct + if (i.lt.inct-1) then + write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i), + & DU2orig(i) + else if (i.eq.inct-1) then + write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i) + else + write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i) + endif + enddo +#endif + ghalf(ind+1)=dmorig(innt) + ghalf(ind+2)=du1orig(innt) + ghalf(ind+3)=dmorig(innt+1) + ind=ind+3 + do i=3,n + ind=ind+i-3 +c write (iout,*) "i",i," ind",ind," indu2",innt+i-2, +c & " indu1",innt+i-1," indm",innt+i + ghalf(ind+1)=du2orig(innt-1+i-2) + ghalf(ind+2)=du1orig(innt-1+i-1) + ghalf(ind+3)=dmorig(innt-1+i) +c write (iout,'(3(a,i2,1x))') "DU2",innt-1+i-2, +c & "DU1",innt-1+i-1,"DM ",innt-1+i + ind=ind+3 + enddo +#ifdef DEBUG + ind=0 + do i=1,n + do j=1,i + ind=ind+1 + inertia(i,j)=ghalf(ind) + inertia(j,i)=ghalf(ind) + enddo + enddo +#endif +#ifdef DEBUG + write (iout,*) "Chain ",ichain," ind",ind," dim",n*(n+1)/2 + write (iout,*) "Five-diagonal inertia matrix, lower triangle" + call matoutr(n,ghalf) +#endif + call gldiag(maxres2_chain,n,n,Ghalf,work,Geigen,Gvec,ierr,iwork) + if (large) then + write (iout,'(//a,i3)') + & "Eigenvectors and eigenvalues of the G matrix chain",ichain + call eigout(n,n,maxres2_chain,maxres2_chain,Gvec,Geigen) + endif +#ifdef DIAGCHECK +c check diagonalization + do i=1,n + do j=1,n + aux=0.0d0 + do k=1,n + do l=1,n + aux=aux+gvec(k,i)*gvec(l,j)*inertia(k,l) + enddo + enddo + if (i.eq.j) then + write (iout,*) i,j,aux,geigen(i) + else + write (iout,*) i,j,aux + endif + enddo + enddo +#endif + xv=0.0d0 + ii=0 + do i=1,n + do k=1,3 + ii=ii+1 + sigv=dsqrt((Rb*t_bath)/geigen(i)) + lowb=-5*sigv + highb=5*sigv + d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb) + EK=EK+0.5d0*geigen(i)*d_t_work_new(ii)**2 +c write (iout,*) "i",i," ii",ii," geigen",geigen(i), +c & " d_t_work_new",d_t_work_new(ii) + enddo + enddo + do k=1,3 + do i=1,n + ind=(i-1)*3+k + d_t_work(ind)=0.0d0 + do j=1,n + d_t_work(ind)=d_t_work(ind) + & +Gvec(i,j)*d_t_work_new((j-1)*3+k) + enddo +c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind) +c call flush(iout) + enddo + enddo +#ifdef DEBUG + aux=0.0d0 + do k=1,3 + do i=1,n + do j=1,n + aux=aux+inertia(i,j)*d_t_work(3*(i-1)+k)*d_t_work(3*(j-1)+k) + enddo + enddo + enddo + Ek3=Ek3+aux/2 +#endif +c Transfer to the d_t vector + innt=chain_border(1,ichain) + inct=chain_border(2,ichain) + ind=0 +c write (iout,*) "ichain",ichain," innt",innt," inct",inct + do i=innt,inct + do j=1,3 + ind=ind+1 + d_t(j,i)=d_t_work(ind) + enddo + 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) + enddo + endif + enddo + enddo + if (large) then + write (iout,*) + write (iout,*) "Random velocities in the Calpha,SC space" + do i=1,nres + write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)') + & restyp(itype(i)),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3) + enddo + endif + call kinetic_CASC(Ek1) +! +! Transform the velocities to virtual-bond space +! +#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 + d_t(j,i)=d_t(j,i+1)-d_t(j,i) + enddo + else + do j=1,3 + d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i) + d_t(j,i)=d_t(j,i+1)-d_t(j,i) + 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=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)) + d_t(:,chain_border1(1,ichain))=0.0d0 + enddo +c write (iout,*) "Adding accelerations" + do ichain=2,nchain +c write (iout,*) "chain",ichain,chain_border1(1,ichain)-1, +c & 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 + 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 +c d_t(j,0)=d_t(j,nnt) +c enddo + do ichain=1,nchain + innt=chain_border(1,ichain) + inct=chain_border(2,ichain) +c write (iout,*) "ichain",ichain," innt",innt," inct",inct +c write (iout,*) "ibeg",ibeg + do j=1,3 + d_t(j,ibeg)=d_t(j,innt) + enddo + ibeg=inct+1 + do i=innt,inct + if (iabs(itype(i).eq.10)) then +c write (iout,*) "i",i,(d_t(j,i),j=1,3),(d_t(j,i+1),j=1,3) + do j=1,3 + d_t(j,i)=d_t(j,i+1)-d_t(j,i) + enddo + else + do j=1,3 + d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i) + d_t(j,i)=d_t(j,i+1)-d_t(j,i) + enddo + end if + enddo + enddo +#endif + if (large) then + write (iout,*) + write (iout,*) + & "Random velocities in the virtual-bond-vector space" + write (iout,'(3hORG,1h(,i5,1h),3f10.5)') 0,(d_t(j,0),j=1,3) + do i=1,nres + write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)') + & restyp(itype(i)),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3) + enddo + write (iout,*) + write (iout,*) "Kinetic energy from inertia matrix eigenvalues", + & Ek + write (iout,*) + & "Kinetic temperatures from inertia matrix eigenvalues", + & 2*Ek/(3*dimen*Rb) +#ifdef DEBUG + write (iout,*) "Kinetic energy from inertia matrix",Ek3 + write (iout,*) "Kinetic temperatures from inertia", + & 2*Ek3/(3*dimen*Rb) +#endif + write (iout,*) "Kinetic energy from velocities in CA-SC space", + & Ek1 + write (iout,*) + & "Kinetic temperatures from velovities in CA-SC space", + & 2*Ek1/(3*dimen*Rb) + call kinetic(Ek1) + write (iout,*) + & "Kinetic energy from virtual-bond-vector velocities",Ek1 + write (iout,*) + & "Kinetic temperature from virtual-bond-vector velocities ", + & 2*Ek1/(dimen3*Rb) + endif +#else xv=0.0d0 ii=0 do i=1,dimen @@ -2052,13 +2560,14 @@ c call kinetic(EK) c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature", c & 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1 c call flush(iout) +#endif return end #ifndef LANG0 c----------------------------------------------------------- subroutine sd_verlet_p_setup c Sets up the parameters of stochastic Verlet algorithm - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' #ifdef MPI include 'mpif.h' @@ -2069,8 +2578,12 @@ c Sets up the parameters of stochastic Verlet algorithm #ifndef LANG0 include 'COMMON.LANGEVIN' #else +#ifdef FIVEDIAG + include 'COMMON.LANGEVIN.lang0.5diag' +#else include 'COMMON.LANGEVIN.lang0' #endif +#endif include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' @@ -2179,14 +2692,12 @@ c c c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables c -#ifndef LANG0 call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat) call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat) call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat) call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat) call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1) call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2) -#endif #ifdef MPI t_sdsetup=t_sdsetup+MPI_Wtime() #else @@ -2229,16 +2740,25 @@ c------------------------------------------------------------- c------------------------------------------------------------- subroutine sd_verlet1 c Applying stochastic velocity Verlet algorithm - step 1 to velocities - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' +#ifdef FIVEDIAG + include 'COMMON.LAGRANGE.5diag' +#else + include 'COMMON.LAGRANGE' +#endif #ifndef LANG0 include 'COMMON.LANGEVIN' #else +#ifdef FIVEDIAG + include 'COMMON.LANGEVIN.lang0.5diag' +#else include 'COMMON.LANGEVIN.lang0' #endif +#endif include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' @@ -2249,6 +2769,7 @@ c Applying stochastic velocity Verlet algorithm - step 1 to velocities double precision stochforcvec(MAXRES6) common /stochcalc/ stochforcvec logical lprn /.false./ + integer i,j,ind,inres c write (iout,*) "dc_old" c do i=0,nres @@ -2309,8 +2830,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) @@ -2332,16 +2853,25 @@ c enddo c-------------------------------------------------------------------------- subroutine sd_verlet2 c Calculating the adjusted velocities for accelerations - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' +#ifdef FIVEDIAG + include 'COMMON.LAGRANGE.5diag' +#else + include 'COMMON.LAGRANGE' +#endif #ifndef LANG0 include 'COMMON.LANGEVIN' #else +#ifdef FIVEDIAG + include 'COMMON.LANGEVIN.lang0.5diag' +#else include 'COMMON.LANGEVIN.lang0' #endif +#endif include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' @@ -2351,6 +2881,7 @@ c Calculating the adjusted velocities for accelerations include 'COMMON.NAMES' double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6) common /stochcalc/ stochforcvec + integer i,j,ind,inres c c Compute the stochastic forces which contribute to velocity change c @@ -2393,7 +2924,7 @@ c----------------------------------------------------------- subroutine sd_verlet_ciccotti_setup c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's c version - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' #ifdef MPI include 'mpif.h' @@ -2401,11 +2932,12 @@ c version include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' -#ifndef LANG0 - include 'COMMON.LANGEVIN' +#ifdef FIVEDIAG + include 'COMMON.LAGRANGE.5diag' #else - include 'COMMON.LANGEVIN.lang0' + include 'COMMON.LAGRANGE' #endif + include 'COMMON.LANGEVIN' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' @@ -2422,6 +2954,7 @@ c version logical lprn /.false./ double precision zero /1.0d-8/, gdt_radius /0.05d0/ double precision ktm + integer i #ifdef MPI tt0 = MPI_Wtime() #else @@ -2492,7 +3025,7 @@ c c------------------------------------------------------------- subroutine sd_verlet1_ciccotti c Applying stochastic velocity Verlet algorithm - step 1 to velocities - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' #ifdef MPI include 'mpif.h' @@ -2500,11 +3033,12 @@ c Applying stochastic velocity Verlet algorithm - step 1 to velocities include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' -#ifndef LANG0 - include 'COMMON.LANGEVIN' +#ifdef FIVEDIAG + include 'COMMON.LAGRANGE.5diag' #else - include 'COMMON.LANGEVIN.lang0' + include 'COMMON.LAGRANGE' #endif + include 'COMMON.LANGEVIN' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' @@ -2515,6 +3049,7 @@ c Applying stochastic velocity Verlet algorithm - step 1 to velocities double precision stochforcvec(MAXRES6) common /stochcalc/ stochforcvec logical lprn /.false./ + integer i,j c write (iout,*) "dc_old" c do i=0,nres @@ -2599,16 +3134,17 @@ c enddo c-------------------------------------------------------------------------- subroutine sd_verlet2_ciccotti c Calculating the adjusted velocities for accelerations - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' -#ifndef LANG0 - include 'COMMON.LANGEVIN' +#ifdef FIVEDIAG + include 'COMMON.LAGRANGE.5diag' #else - include 'COMMON.LANGEVIN.lang0' + include 'COMMON.LAGRANGE' #endif + include 'COMMON.LANGEVIN' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' @@ -2618,6 +3154,7 @@ c Calculating the adjusted velocities for accelerations include 'COMMON.NAMES' double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6) common /stochcalc/ stochforcvec + integer i,j c c Compute the stochastic forces which contribute to velocity change c