X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2FMD.F90;h=87f739439c34f9cbe7cab8d2783e3b1dbc3a70da;hb=91b9f78d94b96277537615722323ebe03cc0a014;hp=9c0b7024f5f674f46159bd174177afde417dbb0e;hpb=cb9719281de8ac538081b7bfca54ff4889a1544d;p=unres4.git diff --git a/source/unres/MD.F90 b/source/unres/MD.F90 index 9c0b702..87f7394 100644 --- a/source/unres/MD.F90 +++ b/source/unres/MD.F90 @@ -554,7 +554,7 @@ ! Calculate energy and forces call zerograd call etotal(potEcomp) - potE=potEcomp(0)-potEcomp(20) + potE=potEcomp(0)-potEcomp(51) call cartgrad totT=totT+d_time totTafm=totT @@ -645,6 +645,122 @@ end subroutine gauss !----------------------------------------------------------------------------- ! kinetic_lesyng.f +#ifdef FIVEDIAG + subroutine kinetic(KE_total) +!c---------------------------------------------------------------- +!c This subroutine calculates the total kinetic energy of the chain +!c----------------------------------------------------------------- +!c 3/5/2020 AL Corrected for multichain systems, no fake peptide groups +!c inside, implemented with five-diagonal inertia matrix + use energy_data + implicit none + real(kind=8):: KE_total,KEt_p,KEt_sc,KEr_p,KEr_sc,mag1,mag2 + integer i,j,k,iti,mnum + real(kind=8),dimension(3) :: incr,v + + KEt_p=0.0d0 + KEt_sc=0.0d0 + KEr_p=0.0D0 + KEr_sc=0.0D0 +!c write (iout,*) "ISC",(isc(itype(i)),i=1,nres) +!c The translational part for peptide virtual bonds + do j=1,3 + incr(j)=d_t(j,0) + enddo + do i=nnt,nct-1 !czy na pewno nct-1?? + mnum=molnum(i) +!c write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3 +!c Skip dummy peptide groups + if (itype(i,mnum).ne.ntyp1_molec(mnum)& + .and. itype(i+1,mnum).ne.ntyp1_molec(mnum)) then + do j=1,3 + v(j)=incr(j)+0.5d0*d_t(j,i) + enddo + if (mnum.eq.5) mp(mnum)=0.0d0 +! if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum) +!c write (iout,*) "Kinetic trp:",i,(v(j),j=1,3) + vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3) + KEt_p=KEt_p+mp(mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3)) + endif + do j=1,3 + incr(j)=incr(j)+d_t(j,i) + enddo + enddo +!c write(iout,*) 'KEt_p', KEt_p +!c The translational part for the side chain virtual bond +!c Only now we can initialize incr with zeros. It must be equal +!c to the velocities of the first Calpha. + do j=1,3 + incr(j)=d_t(j,0) + enddo + do i=nnt,nct + mnum=molnum(i) + iti=iabs(itype(i,mnum)) + if (mnum.eq.5) iti=itype(i,mnum) + if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)& + .or.mnum.ge.3) then + do j=1,3 + v(j)=incr(j) + enddo + else + do j=1,3 + v(j)=incr(j)+d_t(j,nres+i) + enddo + endif +! if (mnum.ne.5) then +! write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3) +! write (iout,*) "i",i," msc",msc(iti,mnum)," v",(v(j),j=1,3) + KEt_sc=KEt_sc+msc(iti,mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3)) + vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3) +! endif + do j=1,3 + incr(j)=incr(j)+d_t(j,i) + enddo + enddo +! goto 111 +! write(iout,*) 'KEt_sc', KEt_sc +! The part due to stretching and rotation of the peptide groups + do i=nnt,nct-1 + mnum=molnum(i) + if (itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.itype(i+1,mnum).ne.ntyp1_molec(mnum)) then + if (mnum.eq.5) Ip(mnum)=Isc(itype(i,mnum),mnum) +! write (iout,*) "i",i +! write (iout,*) "i",i," mag1",mag1," mag2",mag2 + do j=1,3 + incr(j)=d_t(j,i) + enddo +!c write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3) + KEr_p=KEr_p+Ip(mnum)*(incr(1)*incr(1)+incr(2)*incr(2) & + +incr(3)*incr(3)) + endif + enddo +!c goto 111 +!c write(iout,*) 'KEr_p', KEr_p +!c The rotational part of the side chain virtual bond + do i=nnt,nct + mnum=molnum(i) + iti=iabs(itype(i,mnum)) + if (itype(i,1).ne.10.and.itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.mnum.lt.3) then + do j=1,3 + incr(j)=d_t(j,nres+i) + enddo +! write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3) + KEr_sc=KEr_sc+Isc(iti,mnum)*(incr(1)*incr(1)+incr(2)*incr(2)+& + incr(3)*incr(3)) + endif + enddo +!c The total kinetic energy + 111 continue +! write(iout,*) ' KEt_p',KEt_p,' KEt_sc',KEt_sc,' KEr_p',KEr_p, & +! ' KEr_sc', KEr_sc + KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc) +!c write (iout,*) "KE_total",KE_total + return + end subroutine kinetic +#else + !----------------------------------------------------------------------------- subroutine kinetic(KE_total) !---------------------------------------------------------------- @@ -663,7 +779,7 @@ ! include 'COMMON.IOUNITS' real(kind=8) :: KE_total,mscab - integer :: i,j,k,iti,mnum + integer :: i,j,k,iti,mnum,term real(kind=8) :: KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),& mag1,mag2,v(3) #ifdef DEBUG @@ -680,13 +796,21 @@ do j=1,3 incr(j)=d_t(j,0) enddo - do i=nnt,nct-1 + term=nct-1 +! if (molnum(nct).gt.3) term=nct + do i=nnt,term mnum=molnum(i) - if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum) -! write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3) + if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum) +! write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3),mp(mnum) + if (mnum.gt.4) then do j=1,3 v(j)=incr(j)+0.5d0*d_t(j,i) - enddo + enddo + else + do j=1,3 + v(j)=incr(j)+0.5d0*d_t(j,i) + enddo + endif vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3) KEt_p=KEt_p+mp(mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3)) do j=1,3 @@ -703,14 +827,14 @@ do i=nnt,nct mnum=molnum(i) iti=iabs(itype(i,mnum)) - if (mnum.eq.5) then - mscab=0.0d0 - else +! if (mnum.ge.4) then +! mscab=0.0d0 +! else mscab=msc(iti,mnum) - endif +! endif ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)& - .or.(mnum.eq.5)) then + .or.(mnum.ge.4)) then do j=1,3 v(j)=incr(j) enddo @@ -720,7 +844,7 @@ enddo endif ! write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3) -! write (iout,*) "i",i," msc",msc(iti)," v",(v(j),j=1,3) +! write (iout,*) "i",i," msc",msc(iti,mnum)," v",(v(j),j=1,3) KEt_sc=KEt_sc+mscab*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3)) vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3) do j=1,3 @@ -751,7 +875,7 @@ iti=iabs(itype(i,mnum)) ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& - .and.(mnum.ne.5)) then + .and.(mnum.lt.4)) then do j=1,3 incr(j)=d_t(j,nres+i) enddo @@ -768,6 +892,110 @@ return end subroutine kinetic !----------------------------------------------------------------------------- +#endif + subroutine kinetic_CASC(KE_total) +!c---------------------------------------------------------------- +!c Compute the kinetic energy of the system using the Calpha-SC +!c coordinate system +!c----------------------------------------------------------------- + implicit none + double precision KE_total + + integer i,j,k,iti,ichain,innt,inct,mnum + double precision KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),& + mag1,mag2,v(3) +#ifdef FIVEDIAG + KEt_p=0.0d0 + KEt_sc=0.0d0 + KEr_p=0.0D0 + KEr_sc=0.0D0 +!c write (iout,*) "ISC",(isc(itype(i)),i=1,nres) +!c The translational part for peptide virtual bonds + do ichain=1,nchain + + innt=chain_border(1,ichain) + inct=chain_border(2,ichain) +!c write (iout,*) "Kinetic_CASC chain",ichain," innt",innt, +!c & " inct",inct + + do i=innt,inct-1 + mnum=molnum(i) + if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum) + if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum) +!c write (iout,*) i,(d_t(j,i),j=1,3),(d_t(j,i+1),j=1,3) + do j=1,3 + v(j)=0.5d0*(d_t(j,i)+d_t(j,i+1)) + enddo +!c write (iout,*) "Kinetic trp i",i," v",(v(j),j=1,3) + KEt_p=KEt_p+mp(mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3)) + enddo +!c write(iout,*) 'KEt_p', KEt_p +!c The translational part for the side chain virtual bond +!c Only now we can initialize incr with zeros. It must be equal +!c to the velocities of the first Calpha. + do i=innt,inct + mnum=molnum(i) + iti=iabs(itype(i,mnum)) + if (itype(i,1).eq.10.or.mnum.ge.3.or. itype(i,mnum).eq.ntyp1_molec(mnum)) then +!c write (iout,*) i,iti,(d_t(j,i),j=1,3) + do j=1,3 + v(j)=d_t(j,i) + enddo + else +!c write (iout,*) i,iti,(d_t(j,nres+i),j=1,3) + do j=1,3 + v(j)=d_t(j,nres+i) + enddo + endif +!c write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3) +!c write (iout,*) "i",i," msc",msc(iti)," v",(v(j),j=1,3) + KEt_sc=KEt_sc+msc(iti,mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3)) + enddo +!c goto 111 +!c write(iout,*) 'KEt_sc', KEt_sc +!c The part due to stretching and rotation of the peptide groups + do i=innt,inct-1 + mnum=molnum(i) + do j=1,3 + incr(j)=d_t(j,i+1)-d_t(j,i) + enddo + if (mnum.eq.5) Ip(mnum)=Isc(itype(i,mnum),mnum) +!c write (iout,*) i,(incr(j),j=1,3) +!c write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3) + KEr_p=KEr_p+Ip(mnum)*(incr(1)*incr(1)+incr(2)*incr(2)& + +incr(3)*incr(3)) + enddo +!c goto 111 +!c write(iout,*) 'KEr_p', KEr_p +!c The rotational part of the side chain virtual bond + do i=innt,inct + mnum=molnum(i) + iti=iabs(itype(i,mnum)) +! if (iti.ne.10.and.mnum.lt.3) then + if (itype(i,1).ne.10.and.mnum.lt.3.and. itype(i,mnum).ne.ntyp1_molec(mnum)) then + do j=1,3 + incr(j)=d_t(j,nres+i)-d_t(j,i) + enddo +!c write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3) + KEr_sc=KEr_sc+Isc(iti,mnum)*(incr(1)*incr(1)+incr(2)*incr(2)+& + incr(3)*incr(3)) + endif + enddo + + enddo ! ichain +!c The total kinetic energy + 111 continue +!c write(iout,*) ' KEt_p',KEt_p,' KEt_sc',KEt_sc,' KEr_p',KEr_p, +!c & ' KEr_sc', KEr_sc + KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc) +!c write (iout,*) "KE_total",KE_tota +#else + write (iout,*) "Need to compile with -DFIVEDIAG to use this sub!" + stop +#endif + return + end subroutine kinetic_CASC + ! MD_A-MTS.F !----------------------------------------------------------------------------- subroutine MD @@ -814,7 +1042,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 @@ -904,6 +1132,7 @@ tt0=tcpu() #endif do itime=1,n_timestep + if (large) print *,itime,ntwe if (ovrtim()) exit if (large.and. mod(itime,ntwe).eq.0) & write (iout,*) "itime",itime @@ -934,7 +1163,9 @@ enddo #endif else if (lang.eq.1 .or. lang.eq.4) then + print *,"before setup_fricmat" call setup_fricmat + print *,"after setup_fricmat" endif write (iout,'(a,i10)') & "Friction matrix reset based on surface area, itime",itime @@ -973,7 +1204,9 @@ call RESPA_step(itime) else ! Variable time step algorithm. + if (large) print *,"before verlet_step" call velverlet_step(itime) + if (large) print *,"after verlet_step" endif else #ifdef BROWN @@ -989,8 +1222,9 @@ itime_mat=itime if (ntwe.ne.0) then if (mod(itime,ntwe).eq.0) then +! call returnbox call statout(itime) - call returnbox +! call returnbox ! call check_ecartint endif #ifdef VOUT @@ -1006,7 +1240,7 @@ enddo do i=nnt,nct mnum=molnum(i) - if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.mnum.ne.5) then + if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum).and.mnum.lt.4) then do j=1,3 ind=ind+1 v_work(ind)=d_t(j,i+nres) @@ -1028,14 +1262,31 @@ endif if (mod(itime,ntwx).eq.0) then call returnbox + call enerprint(potEcomp) + write (tytul,'("time",f8.2)') totT if(mdpdb) then call hairpin(.true.,nharp,iharp) call secondary2(.true.) call pdbout(potE,tytul,ipdb) +! call enerprint(potEcomp) else call cartout(totT) endif + if (fodson) then + write(iout,*) "starting fodstep" + call fodstep(nfodstep) + write(iout,*) "after fodstep" + call statout(itime) + if(mdpdb) then + call hairpin(.true.,nharp,iharp) + call secondary2(.true.) + call pdbout(potE,tytul,ipdb) + else + call cartout(totT) + endif + endif + endif if (rstcount.eq.1000.or.itime.eq.n_timestep) then open(irest2,file=rest2name,status='unknown') @@ -1130,6 +1381,7 @@ icount_scale=0 if (lang.eq.1) then call sddir_precalc + if (large) print *,"after sddir_precalc" else if (lang.eq.2 .or. lang.eq.3) then #ifndef LANG0 call stochastic_force(stochforcvec) @@ -1176,7 +1428,7 @@ call chainbuild_cart if (rattle) call rattle1 if (ntwe.ne.0) then - if (large.and. mod(itime,ntwe).eq.0) then + if (large) then !.and. mod(itime,ntwe).eq.0) then write (iout,*) "Cartesian and internal coordinates: step 1" call cartprint call intout @@ -1232,7 +1484,7 @@ t_etotal=t_etotal+tcpu()-tt0 #endif #endif - potE=potEcomp(0)-potEcomp(20) + potE=potEcomp(0)-potEcomp(51) call cartgrad ! Get the new accelerations call lagrangian @@ -1747,6 +1999,7 @@ t_elong=t_elong+tcpu()-tt0 #endif #endif + potE=potEcomp(0)-potEcomp(51) call cartgrad call lagrangian #ifdef MPI @@ -1781,7 +2034,7 @@ do i=0,n_ene potEcomp(i)=energia_short(i)+energia_long(i) enddo - potE=potEcomp(0)-potEcomp(20) + potE=potEcomp(0)-potEcomp(51) ! potE=energia_short(0)+energia_long(0) totT=totT+d_time totTafm=totT @@ -1836,7 +2089,7 @@ mnum=molnum(i) ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& - .and.(mnum.ne.5)) then + .and.(mnum.lt.4)) then inres=i+nres do j=1,3 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time @@ -1891,7 +2144,7 @@ mnum=molnum(i) ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& - .and.(mnum.ne.5)) then + .and.(mnum.lt.4)) then inres=i+nres do j=1,3 adt=d_a_old(j,inres)*d_time @@ -1942,7 +2195,7 @@ ! iti=iabs(itype(i,mnum)) ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& - .and.(mnum.ne.5)) then + .and.(mnum.lt.4)) then inres=i+nres do j=1,3 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time @@ -1980,12 +2233,15 @@ !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres !el common /stochcalc/ stochforcvec real(kind=8) :: time00 + integer :: i ! ! Compute friction and stochastic forces ! #ifdef MPI time00=MPI_Wtime() + if (large) print *,"before friction_force" call friction_force + if (large) print *,"after friction_force" time_fric=time_fric+MPI_Wtime()-time00 time00=MPI_Wtime() call stochastic_force(stochforcvec) @@ -1995,8 +2251,23 @@ ! Compute the acceleration due to friction forces (d_af_work) and stochastic ! forces (d_as_work) ! +! call ginv_mult(fric_work, d_af_work) +! call ginv_mult(stochforcvec, d_as_work) +#ifdef FIVEDIAG + call fivediaginv_mult(dimen,fric_work, d_af_work) + call fivediaginv_mult(dimen,stochforcvec, d_as_work) + if (large) then + write(iout,*),"dimen",dimen + do i=1,dimen + write(iout,*) "fricwork",fric_work(i), d_af_work(i) + write(iout,*) "stochforcevec", stochforcvec(i), d_as_work(i) + enddo + endif +#else call ginv_mult(fric_work, d_af_work) call ginv_mult(stochforcvec, d_as_work) +#endif + return end subroutine sddir_precalc !----------------------------------------------------------------------------- @@ -2033,6 +2304,7 @@ do j=1,3 adt=(d_a_old(j,0)+d_af_work(j))*d_time adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time +! write(iout,*) i,"adt",adt,"ads",adt2,d_a_old(j,0),d_af_work(j),d_time dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt d_t(j,0)=d_t_old(j,0)+adt @@ -2042,6 +2314,7 @@ 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 +! write(iout,*) i,"adt",adt,"ads",adt2,d_a_old(j,i),d_af_work(ind+j) dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt d_t(j,i)=d_t_old(j,i)+adt @@ -2053,11 +2326,12 @@ ! iti=iabs(itype(i,mnum)) ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& - .and.(mnum.ne.5)) then + .and.(mnum.lt.4)) then inres=i+nres 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,*) i,"adt",adt,"ads",adt2,d_a_old(j,inres),d_af_work(ind+j) 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 @@ -2065,6 +2339,7 @@ ind=ind+3 endif enddo + return end subroutine sddir_verlet1 !----------------------------------------------------------------------------- @@ -2103,7 +2378,11 @@ ! Compute the acceleration due to friction forces (d_af_work) and stochastic ! forces (d_as_work) ! +#ifdef FIVEDIAG + call fivediaginv_mult(6*nres,stochforcvec, d_as_work1) +#else call ginv_mult(stochforcvec, d_as_work1) +#endif ! ! Update velocities @@ -2125,7 +2404,7 @@ ! iti=iabs(itype(i,mnum)) ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& - .and.(mnum.ne.5)) then + .and.(mnum.lt.4)) then inres=i+nres do j=1,3 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres) & @@ -2199,7 +2478,7 @@ ! iti=iabs(itype(i,mnum)) ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& - .and.(mnum.ne.5)) then + .and.(mnum.lt.4)) then do j=1,3 ! accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres) accel_old(j)=accel_old(j)+d_a_old(j,i+nres) @@ -2257,7 +2536,7 @@ endif ! Side chains if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.& - molnum(i).ne.5) then + molnum(i).lt.4) then do j=1,3 epdriftij= & dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i)) @@ -2310,7 +2589,7 @@ ! iti=iabs(itype(i,mnum)) ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& - .and.(mnum.ne.5)) then + .and.(mnum.lt.4)) then inres=i+nres do j=1,3 d_t(j,inres)=fact*d_t(j,inres) @@ -2331,6 +2610,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 +2618,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 +2639,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) @@ -2507,7 +2790,7 @@ if(me.eq.king.or..not.out1file)then write (iout,*) "Initial velocities" do i=0,nres - write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),& + write (iout,'(i6,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),& (d_t(j,i+nres),j=1,3) enddo ! Zeroing the total angular momentum of the system @@ -2527,8 +2810,14 @@ write (iout,*) "vcm right after adjustment:" write (iout,*) (vcm(j),j=1,3) endif - if (.not.rest) then + + + ! 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' @@ -2552,9 +2841,13 @@ call minimize(etot,varia,iretcode,nfun) call var_to_geom(nvar,varia) endif + write(iout,*) "just before minimin" + call cartprint if(me.eq.king.or..not.out1file) & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun endif + write(iout,*) "just after minimin" + call cartprint if(iranconf.ne.0) then !c 8/22/17 AL Loop to produce a low-energy random conformation DO iranmin=1,40 @@ -2583,7 +2876,8 @@ endif if(me.eq.king.or..not.out1file) & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun - + write(iout,*) "just after minimin" + call cartprint if (isnan(etot) .or. etot.gt.4.0d6) then write (iout,*) "Energy too large",etot, & " trying another random conformation" @@ -2633,15 +2927,148 @@ 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 call kinetic(EK) + write(iout,*) "just after kinetic" + call cartprint if (tbf) then call verlet_bath endif kinetic_T=2.0d0/(dimen3*Rb)*EK if(me.eq.king.or..not.out1file)then + write(iout,*) "just after verlet_bath" call cartprint call intout endif @@ -2663,7 +3090,9 @@ #endif potE=potEcomp(0) call cartgrad + write(iout,*) "before lagrangian" call lagrangian + write(iout,*) "before max_accel" call max_accel if (amax*d_time .gt. dvmax) then d_time=d_time*dvmax/amax @@ -2676,7 +3105,7 @@ call enerprint(potEcomp) ! write(iout,*) (potEcomp(i),i=0,n_ene) endif - potE=potEcomp(0)-potEcomp(20) + potE=potEcomp(0)-potEcomp(51) totE=EK+potE itime=0 itime_mat=itime @@ -2700,7 +3129,7 @@ write (iout,*) "Initial velocities" write (iout,"(13x,' backbone ',23x,' side chain')") do i=0,nres - write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),& + write (iout,'(i6,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),& (d_t(j,i+nres),j=1,3) enddo write (iout,*) "Initial accelerations" @@ -2810,318 +3239,267 @@ ! include 'COMMON.NAMES' ! include 'COMMON.TIME1' real(kind=8) :: xv,sigv,lowb,highb ,Ek1 +#ifdef FIVEDIAG + integer ichain,n,innt,inct,ibeg,ierr + real(kind=8) ,allocatable, dimension(:):: work + integer,allocatable,dimension(:) :: iwork +! 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 +#endif !#define DEBUG #ifdef FIVEDIAG - real(kind=8) ,allocatable, dimension(:) :: DDU1,DDU2,DL2,DL1,xsolv,DML,rs - real(kind=8) :: sumx + real(kind=8) ,allocatable, dimension(:) :: xsolv,DML,rs + real(kind=8) :: sumx,Ek2,Ek3,aux #ifdef DEBUG real(kind=8) ,allocatable, dimension(:) :: rsold - real (kind=8),allocatable,dimension(:,:) :: matold + real (kind=8),allocatable,dimension(:,:) :: matold,inertia integer :: iti #endif #endif - integer :: i,j,ii,k,ind,mark,imark,mnum + integer :: i,j,ii,k,mark,imark,mnum,nres2 + integer(kind=8) :: ind ! Generate random velocities from Gaussian distribution of mean 0 and std of KT/m +!#undef DEBUG ! First generate velocities in the eigenspace of the G matrix ! write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3 ! call flush(iout) - xv=0.0d0 - ii=0 - do i=1,dimen - 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) +#ifdef FIVEDIAG + if(.not.allocated(work)) then + allocate(work(48*nres)) + allocate(iwork(6*nres)) + endif + print *,"IN RANDOM VEL" + nres2=2*nres +! print *,size(ghalf) +#undef DEBUG #ifdef DEBUG - write (iout,*) "i",i," ii",ii," geigen",geigen(i),& - " d_t_work_new",d_t_work_new(ii) + write (iout,*) "Random_vel, fivediag" + flush(iout) + allocate(inertia(2*nres,2*nres)) #endif - enddo - enddo + d_t=0.0d0 + Ek2=0.0d0 + EK=0.0d0 + Ek3=0.0d0 #ifdef DEBUG -! diagnostics - Ek1=0.0d0 - ii=0 - do i=1,dimen - do k=1,3 - ii=ii+1 - Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2 - enddo - enddo - write (iout,*) "Ek from eigenvectors",Ek1 - write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb) -! end diagnostics + write(iout,*), nchain #endif -#ifdef FIVEDIAG -! Transform velocities to UNRES coordinate space - allocate (DL1(2*nres)) - allocate (DDU1(2*nres)) - allocate (DL2(2*nres)) - allocate (DDU2(2*nres)) - allocate (xsolv(2*nres)) - allocate (DML(2*nres)) - allocate (rs(2*nres)) + do ichain=1,nchain + ind=0 + if(.not.allocated(ghalf)) print *,"COCO" + if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) + ghalf=0.0d0 + n=dimen_chain(ichain) + innt=iposd_chain(ichain) + inct=innt+n-1 #ifdef DEBUG - allocate (rsold(2*nres)) - allocate (matold(2*nres,2*nres)) - matold=0 - matold(1,1)=DMorig(1) - matold(1,2)=DU1orig(1) - matold(1,3)=DU2orig(1) - write (*,*) DMorig(1),DU1orig(1),DU2orig(1) - - do i=2,dimen-1 - do j=1,dimen - if (i.eq.j) then - matold(i,j)=DMorig(i) - matold(i,j-1)=DU1orig(i-1) - if (j.ge.3) then - matold(i,j-2)=DU2orig(i-2) - endif - - endif - - enddo - do j=1,dimen-1 - if (i.eq.j) then - matold(i,j+1)=DU1orig(i) - - end if - enddo - do j=1,dimen-2 - if(i.eq.j) then - matold(i,j+2)=DU2orig(i) - endif - enddo - enddo - matold(dimen,dimen)=DMorig(dimen) - matold(dimen,dimen-1)=DU1orig(dimen-1) - matold(dimen,dimen-2)=DU2orig(dimen-2) - write(iout,*) "old gmatrix" - call matout(dimen,dimen,2*nres,2*nres,matold) + 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 - d_t_work=0.0d0 - do i=1,dimen -! Find the ith eigenvector of the pentadiagonal inertiq matrix - do imark=dimen,1,-1 - do j=1,imark-1 - DML(j)=DMorig(j)-geigen(i) - enddo - do j=imark+1,dimen - DML(j-1)=DMorig(j)-geigen(i) - enddo - do j=1,imark-2 - DDU1(j)=DU1orig(j) - enddo - DDU1(imark-1)=DU2orig(imark-1) - do j=imark+1,dimen-1 - DDU1(j-1)=DU1orig(j) - enddo - do j=1,imark-3 - DDU2(j)=DU2orig(j) - enddo - DDU2(imark-2)=0.0d0 - DDU2(imark-1)=0.0d0 - do j=imark,dimen-3 - DDU2(j)=DU2orig(j+1) - enddo - do j=1,dimen-3 - DL2(j+2)=DDU2(j) - enddo - do j=1,dimen-2 - DL1(j+1)=DDU1(j) - enddo -#ifdef DEBUG - write (iout,*) "DL2,DL1,DML,DDU1,DDU2" - write(iout,'(10f10.5)') (DL2(k),k=3,dimen-1) - write(iout,'(10f10.5)') (DL1(k),k=2,dimen-1) - write(iout,'(10f10.5)') (DML(k),k=1,dimen-1) - write(iout,'(10f10.5)') (DDU1(k),k=1,dimen-2) - write(iout,'(10f10.5)') (DDU2(k),k=1,dimen-3) -#endif - rs=0.0d0 - if (imark.gt.2) rs(imark-2)=-DU2orig(imark-2) - if (imark.gt.1) rs(imark-1)=-DU1orig(imark-1) - if (imark.lt.dimen) rs(imark)=-DU1orig(imark) - if (imark.lt.dimen-1) rs(imark+1)=-DU2orig(imark) -#ifdef DEBUG - rsold=rs -#endif -! write (iout,*) "Vector rs" -! do j=1,dimen-1 -! write (iout,*) j,rs(j) -! enddo - xsolv=0.0d0 - call FDIAG(dimen-1,DL2,DL1,DML,DDU1,DDU2,rs,xsolv,mark) - - if (mark.eq.1) then - -#ifdef DEBUG -! Check solution - do j=1,imark-1 - sumx=-geigen(i)*xsolv(j) - do k=1,imark-1 - sumx=sumx+matold(j,k)*xsolv(k) - enddo - do k=imark+1,dimen - sumx=sumx+matold(j,k)*xsolv(k-1) - enddo - write(iout,'(i5,3f10.5)') j,sumx,rsold(j),sumx-rsold(j) - enddo - do j=imark+1,dimen - sumx=-geigen(i)*xsolv(j-1) - do k=1,imark-1 - sumx=sumx+matold(j,k)*xsolv(k) - enddo - do k=imark+1,dimen - sumx=sumx+matold(j,k)*xsolv(k-1) - enddo - write(iout,'(i5,3f10.5)') j-1,sumx,rsold(j-1),sumx-rsold(j-1) - enddo -! end check - write (iout,*)& - "Solution of equations system",i," for eigenvalue",geigen(i) - do j=1,dimen-1 - write(iout,'(i5,f10.5)') j,xsolv(j) - enddo -#endif - do j=dimen-1,imark,-1 - xsolv(j+1)=xsolv(j) - enddo - xsolv(imark)=1.0d0 + 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 + write (iout,*) "i",i," ind",ind," indu2",innt+i-2,& + " 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 - write (iout,*) "Un-normalized eigenvector",i," for eigenvalue",geigen(i) - do j=1,dimen - write(iout,'(i5,f10.5)') j,xsolv(j) - enddo + 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 -! Normalize ith eigenvector - sumx=0.0d0 - do j=1,dimen - sumx=sumx+xsolv(j)**2 - enddo - sumx=dsqrt(sumx) - do j=1,dimen - xsolv(j)=xsolv(j)/sumx - enddo #ifdef DEBUG - write (iout,*) "Eigenvector",i," for eigenvalue",geigen(i) - do j=1,dimen - write(iout,'(i5,3f10.5)') j,xsolv(j) - enddo + write (iout,*) "Chain ",ichain," ind",ind," dim",n*(n+1)/2 + write (iout,*) "Five-diagonal inertia matrix, lower triangle" +! call matoutr(n,ghalf) #endif -! All done at this point for eigenvector i, exit loop - exit - - endif ! mark.eq.1 - - enddo ! imark - - if (mark.ne.1) then - write (iout,*) "Unable to find eigenvector",i - endif - -! write (iout,*) "i=",i - do k=1,3 -! write (iout,*) "k=",k - do j=1,dimen - ind=(j-1)*3+k -! write(iout,*) "ind",ind," ind1",3*(i-1)+k - d_t_work(ind)=d_t_work(ind) & - +xsolv(j)*d_t_work_new(3*(i-1)+k) - enddo - enddo - enddo ! i (loop over the eigenvectors) - -#ifdef DEBUG - write (iout,*) "d_t_work" - do i=1,3*dimen - write (iout,"(i5,f10.5)") i,d_t_work(i) - enddo - Ek1=0.0d0 - ii=0 - do i=nnt,nct -! if (itype(i,1).eq.10) then - mnum=molnum(i) - if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum) - iti=iabs(itype(i,mnum)) -! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then - if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)& - .or.(mnum.eq.5)) then - j=ii+3 - else - j=ii+6 + call gldiag(nres*2,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,nres*2,nres*2,Gvec,Geigen) endif - if (i.lt.nct) then +#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 -! write (iout,*) "k",k," ii+k",ii+k," ii+j+k",ii+j+k,"EK1",Ek1 - Ek1=Ek1+0.5d0*mp(mnum)*((d_t_work(ii+j+k)+d_t_work_new(ii+k))/2)**2+& - 0.5d0*0.25d0*IP(mnum)*(d_t_work(ii+j+k)-d_t_work_new(ii+k))**2 + 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 - endif - if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& - .and.(mnum.ne.5)) ii=ii+3 - write (iout,*) "i",i," itype",itype(i,mnum)," mass",msc(itype(i,mnum),mnum) - write (iout,*) "ii",ii + enddo do k=1,3 - ii=ii+1 - write (iout,*) "k",k," ii",ii,"EK1",EK1 - if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& - .and.(mnum.ne.5))& - Ek1=Ek1+0.5d0*Isc(iabs(itype(i,mnum)),mnum)*(d_t_work(ii)-d_t_work(ii-3))**2 - Ek1=Ek1+0.5d0*msc(iabs(itype(i,mnum)),mnum)*d_t_work(ii)**2 + 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 - write (iout,*) "i",i," ii",ii - enddo - write (iout,*) "Ek from d_t_work",Ek1 - write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb) -#endif - if(allocated(DDU1)) deallocate(DDU1) - if(allocated(DDU2)) deallocate(DDU2) - if(allocated(DL2)) deallocate(DL2) - if(allocated(DL1)) deallocate(DL1) - if(allocated(xsolv)) deallocate(xsolv) - if(allocated(DML)) deallocate(DML) - if(allocated(rs)) deallocate(rs) #ifdef DEBUG - if(allocated(matold)) deallocate(matold) - if(allocated(rsold)) deallocate(rsold) -#endif - ind=1 - do j=nnt,nct + aux=0.0d0 do k=1,3 - d_t(k,j)=d_t_work(ind) - ind=ind+1 + 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 - mnum=molnum(j) - if (itype(j,1).ne.10 .and. itype(j,mnum).ne.ntyp1_molec(mnum)& - .and.(mnum.ne.5)) then - do k=1,3 - d_t(k,j+nres)=d_t_work(ind) + 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 - end if + mnum=molnum(i) + if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum).and.mnum.le.2) then + do j=1,3 + ind=ind+1 + d_t(j,i+nres)=d_t_work(ind) + enddo + endif + enddo enddo -#ifdef DEBUG - write (iout,*) "Random velocities in the Calpha,SC space" - do i=nnt,nct - write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3) + if (large) then + write (iout,*) + write (iout,*) "Random velocities in the Calpha,SC space" + do i=1,nres + mnum=molnum(i) + write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')& + restyp(itype(i,mnum),mnum),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 + mnum=molnum(i) + if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum).or.mnum.ge.3) 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 - do i=nnt,nct - write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3) + 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 -#endif +#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,0)=d_t(j,nnt) + d_t(j,ibeg)=d_t(j,innt) enddo - do i=nnt,nct -! if (itype(i,1).eq.10) then - mnum=molnum(i) - if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)& - .or.(mnum.eq.5)) then + ibeg=inct+1 + do i=innt,inct + mnum=molnum(i) + if (iabs(itype(i,1).eq.10).or.mnum.ge.3) 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 @@ -3132,19 +3510,71 @@ 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,mnum),mnum),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,*)"Random velocities in the virtual-bond-vector space" - do i=nnt,nct-1 - write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3) + 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 + 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) +#ifdef DEBUG + write (iout,*) "i",i," ii",ii," geigen",geigen(i),& + " d_t_work_new",d_t_work_new(ii) +#endif + enddo enddo - do i=nnt,nct - write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3) +#ifdef DEBUG +! diagnostics + Ek1=0.0d0 + ii=0 + do i=1,dimen + do k=1,3 + ii=ii+1 + Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2 + enddo enddo - call kinetic(Ek1) - write (iout,*) "Ek from d_t_work",Ek1 + write (iout,*) "Ek from eigenvectors",Ek1 write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb) +! end diagnostics #endif -#else + do k=0,2 do i=1,dimen ind=(i-1)*3+k+1 @@ -3172,7 +3602,7 @@ mnum=molnum(i) ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& - .and.(mnum.ne.5)) then + .and.(mnum.lt.4)) then do j=1,3 ind=ind+1 d_t(j,i+nres)=d_t_work(ind) @@ -3185,6 +3615,7 @@ ! 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1 ! call flush(iout) ! write(iout,*) "end init MD" +#undef DEBUG return end subroutine random_vel !----------------------------------------------------------------------------- @@ -3416,7 +3847,7 @@ mnum=molnum(i) ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& - .and.(mnum.ne.5)) then + .and.(mnum.lt.4)) then do j=1,3 dc_work(ind+j)=dc_old(j,i+nres) d_t_work(ind+j)=d_t_old(j,i+nres) @@ -3467,7 +3898,7 @@ mnum=molnum(i) ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& - .and.(mnum.ne.5)) then + .and.(mnum.lt.4)) then inres=i+nres do j=1,3 dc(j,inres)=dc_work(ind+j) @@ -3536,7 +3967,7 @@ mnum=molnum(i) ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& - .and.(mnum.ne.5)) then + .and.(mnum.lt.4)) then inres=i+nres do j=1,3 d_t(j,inres)=d_t_work(ind+j) @@ -3717,8 +4148,8 @@ " vrand_mat2" do i=1,dimen do j=1,dimen - write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),& - vfric_mat(i,j),afric_mat(i,j),& + write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),& + vfric_mat(i,j),afric_mat(i,j),& prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j) enddo enddo @@ -3752,7 +4183,7 @@ mnum=molnum(i) ! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& - .and.(mnum.ne.5)) then + .and.(mnum.lt.4)) then inres=i+nres do j=1,3 dc(j,inres)=dc_work(ind+j) @@ -3806,35 +4237,371 @@ enddo d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2 enddo -#endif +#endif + do j=1,3 + d_t(j,0)=d_t_work(j) + enddo + ind=3 + do i=nnt,nct-1 + do j=1,3 + d_t(j,i)=d_t_work(ind+j) + enddo + ind=ind+3 + enddo + do i=nnt,nct + mnum=molnum(i) + if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.(mnum.lt.4)) +! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then + inres=i+nres + do j=1,3 + d_t(j,inres)=d_t_work(ind+j) + enddo + ind=ind+3 + endif + enddo + return + end subroutine sd_verlet2_ciccotti +#endif +!----------------------------------------------------------------------------- +! moments.f +!----------------------------------------------------------------------------- +#ifdef FIVEDIAG + subroutine inertia_tensor + use comm_gucio + use energy_data + real(kind=8) Im(3,3),Imcp(3,3),pr(3),M_SC,& + eigvec(3,3),Id(3,3),eigval(3),L(3),vp(3),vrot(3),& + vpp(3,0:MAXRES),vs_p(3),pr1(3,3),& + pr2(3,3),pp(3),incr(3),v(3),mag,mag2,M_PEP + integer iti,inres,i,j,k,mnum,mnum1 + do i=1,3 + do j=1,3 + Im(i,j)=0.0d0 + pr1(i,j)=0.0d0 + pr2(i,j)=0.0d0 + enddo + L(i)=0.0d0 + cm(i)=0.0d0 + vrot(i)=0.0d0 + enddo + M_PEP=0.0d0 + +!c caulating the center of the mass of the protein + do i=nnt,nct-1 + mnum=molnum(i) + mnum1=molnum(i+1) + if (itype(i,mnum).eq.ntyp1_molec(mnum)& + .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle +! if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum) + if (mnum.ge.5) mp(mnum)=0.0d0 + M_PEP=M_PEP+mp(mnum) + + do j=1,3 + cm(j)=cm(j)+(c(j,i)+0.5d0*dc(j,i))*mp(mnum) + enddo + enddo +! do j=1,3 +! cm(j)=mp*cm(j) +! enddo + M_SC=0.0d0 + do i=nnt,nct + mnum=molnum(i) + mnum1=molnum(i+1) + iti=iabs(itype(i,mnum)) + if (iti.eq.ntyp1_molec(mnum)) cycle + M_SC=M_SC+msc(iabs(iti),mnum) + inres=i+nres + do j=1,3 + cm(j)=cm(j)+msc(iabs(iti),mnum)*c(j,inres) + enddo + enddo + do j=1,3 + cm(j)=cm(j)/(M_SC+M_PEP) + enddo + do i=nnt,nct-1 + mnum=molnum(i) + mnum1=molnum(i+1) +! if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum) + if (mnum.ge.5) mp(mnum)=0.0d0 + if (itype(i,mnum).eq.ntyp1_molec(mnum)& + .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle + do j=1,3 + pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j) + enddo + Im(1,1)=Im(1,1)+mp(mnum)*(pr(2)*pr(2)+pr(3)*pr(3)) + Im(1,2)=Im(1,2)-mp(mnum)*pr(1)*pr(2) + Im(1,3)=Im(1,3)-mp(mnum)*pr(1)*pr(3) + Im(2,3)=Im(2,3)-mp(mnum)*pr(2)*pr(3) + Im(2,2)=Im(2,2)+mp(mnum)*(pr(3)*pr(3)+pr(1)*pr(1)) + Im(3,3)=Im(3,3)+mp(mnum)*(pr(1)*pr(1)+pr(2)*pr(2)) + enddo + + do i=nnt,nct + mnum=molnum(i) + iti=iabs(itype(i,mnum)) + if (iti.eq.ntyp1_molec(mnum)) cycle + inres=i+nres + do j=1,3 + pr(j)=c(j,inres)-cm(j) + enddo + Im(1,1)=Im(1,1)+msc(iabs(iti),mnum)*(pr(2)*pr(2)+pr(3)*pr(3)) + Im(1,2)=Im(1,2)-msc(iabs(iti),mnum)*pr(1)*pr(2) + Im(1,3)=Im(1,3)-msc(iabs(iti),mnum)*pr(1)*pr(3) + Im(2,3)=Im(2,3)-msc(iabs(iti),mnum)*pr(2)*pr(3) + Im(2,2)=Im(2,2)+msc(iabs(iti),mnum)*(pr(3)*pr(3)+pr(1)*pr(1)) + Im(3,3)=Im(3,3)+msc(iabs(iti),mnum)*(pr(1)*pr(1)+pr(2)*pr(2)) + enddo + do i=nnt,nct-1 + mnum=molnum(i) + mnum1=molnum(i+1) + if (itype(i,mnum).eq.ntyp1_molec(mnum)& + .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle + Im(1,1)=Im(1,1)+Ip(mnum)*(1-dc_norm(1,i)*dc_norm(1,i))*& + vbld(i+1)*vbld(i+1)*0.25d0 + Im(1,2)=Im(1,2)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(2,i))*& + vbld(i+1)*vbld(i+1)*0.25d0 + Im(1,3)=Im(1,3)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(3,i))*& + vbld(i+1)*vbld(i+1)*0.25d0 + Im(2,3)=Im(2,3)+Ip(mnum)*(-dc_norm(2,i)*dc_norm(3,i))*& + vbld(i+1)*vbld(i+1)*0.25d0 + Im(2,2)=Im(2,2)+Ip(mnum)*(1-dc_norm(2,i)*dc_norm(2,i))*& + vbld(i+1)*vbld(i+1)*0.25d0 + Im(3,3)=Im(3,3)+Ip(mnum)*(1-dc_norm(3,i)*dc_norm(3,i))*& + vbld(i+1)*vbld(i+1)*0.25d0 + enddo + do i=nnt,nct + mnum=molnum(i) + mnum1=molnum(i+1) + iti=iabs(itype(i,mnum)) + if (iti.ne.10 .and. iti.ne.ntyp1_molec(mnum).and.mnum.le.2) then + inres=i+nres + Im(1,1)=Im(1,1)+Isc(iti,mnum)*(1-dc_norm(1,inres)*& + dc_norm(1,inres))*vbld(inres)*vbld(inres) + Im(1,2)=Im(1,2)-Isc(iti,mnum)*(dc_norm(1,inres)*& + dc_norm(2,inres))*vbld(inres)*vbld(inres) + Im(1,3)=Im(1,3)-Isc(iti,mnum)*(dc_norm(1,inres)*& + dc_norm(3,inres))*vbld(inres)*vbld(inres) + Im(2,3)=Im(2,3)-Isc(iti,mnum)*(dc_norm(2,inres)*& + dc_norm(3,inres))*vbld(inres)*vbld(inres) + Im(2,2)=Im(2,2)+Isc(iti,mnum)*(1-dc_norm(2,inres)*& + dc_norm(2,inres))*vbld(inres)*vbld(inres) + Im(3,3)=Im(3,3)+Isc(iti,mnum)*(1-dc_norm(3,inres)*& + dc_norm(3,inres))*vbld(inres)*vbld(inres) + endif + enddo + + call angmom(cm,L) + Im(2,1)=Im(1,2) + Im(3,1)=Im(1,3) + Im(3,2)=Im(2,3) + +!c Copng the Im matrix for the djacob subroutine + do i=1,3 + do j=1,3 + Imcp(i,j)=Im(i,j) + Id(i,j)=0.0d0 + enddo + enddo +!c Finding the eigenvectors and eignvalues of the inertia tensor + call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval) + do i=1,3 + if (dabs(eigval(i)).gt.1.0d-15) then + Id(i,i)=1.0d0/eigval(i) + else + Id(i,i)=0.0d0 + endif + enddo + do i=1,3 + do j=1,3 + Imcp(i,j)=eigvec(j,i) + enddo + enddo + do i=1,3 + do j=1,3 + do k=1,3 + pr1(i,j)=pr1(i,j)+Id(i,k)*Imcp(k,j) + enddo + enddo + enddo + do i=1,3 + do j=1,3 + do k=1,3 + pr2(i,j)=pr2(i,j)+eigvec(i,k)*pr1(k,j) + enddo + enddo + enddo +!c Calculating the total rotational velocity of the molecule + do i=1,3 + do j=1,3 + vrot(i)=vrot(i)+pr2(i,j)*L(j) + enddo + enddo + do i=nnt,nct-1 + mnum=molnum(i) + mnum1=molnum(i+1) +! write (iout,*) itype(i+1,mnum1),itype(i,mnum) + if (itype(i+1,mnum1).ne.ntyp1_molec(mnum1) & + .and. itype(i,mnum).eq.ntyp1_molec(mnum) .or.& + itype(i,mnum).ne.ntyp1_molec(mnum) & + .and. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle + call vecpr(vrot(1),dc(1,i),vp) + do j=1,3 + d_t(j,i)=d_t(j,i)-vp(j) + enddo + enddo + do i=nnt,nct + mnum=molnum(i) + if(itype(i,1).ne.10 .and.itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.mnum.le.2) then + inres=i+nres + call vecpr(vrot(1),dc(1,inres),vp) + do j=1,3 + d_t(j,inres)=d_t(j,inres)-vp(j) + enddo + endif + enddo + call angmom(cm,L) + return + end subroutine inertia_tensor +!------------------------------------------------------------ + subroutine angmom(cm,L) + implicit none + double precision L(3),cm(3),pr(3),vp(3),vrot(3),incr(3),v(3),& + pp(3),mscab + integer iti,inres,i,j,mnum,mnum1 +!c Calculate the angular momentum + do j=1,3 + L(j)=0.0d0 + enddo + do j=1,3 + incr(j)=d_t(j,0) + enddo + do i=nnt,nct-1 + mnum=molnum(i) + mnum1=molnum(i+1) +! if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum) + if (mnum.ge.5) mp(mnum)=0.0d0 + if (itype(i,mnum).eq.ntyp1_molec(mnum)& + .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle + do j=1,3 + pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j) + enddo + do j=1,3 + v(j)=incr(j)+0.5d0*d_t(j,i) + enddo + do j=1,3 + incr(j)=incr(j)+d_t(j,i) + enddo + call vecpr(pr(1),v(1),vp) + do j=1,3 + L(j)=L(j)+mp(mnum)*vp(j) + enddo + do j=1,3 + pr(j)=0.5d0*dc(j,i) + pp(j)=0.5d0*d_t(j,i) + enddo + call vecpr(pr(1),pp(1),vp) + do j=1,3 + L(j)=L(j)+Ip(mnum)*vp(j) + enddo + enddo do j=1,3 - d_t(j,0)=d_t_work(j) + incr(j)=d_t(j,0) enddo - ind=3 - do i=nnt,nct-1 + do i=nnt,nct + mnum=molnum(i) + iti=iabs(itype(i,mnum)) + inres=i+nres do j=1,3 - d_t(j,i)=d_t_work(ind+j) + pr(j)=c(j,inres)-cm(j) enddo - ind=ind+3 - enddo - do i=nnt,nct - mnum=molnum(i) - if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& - .and.(mnum.ne.5)) -! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then - inres=i+nres + if (itype(i,1).ne.10 .and.itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.mnum.le.2) then do j=1,3 - d_t(j,inres)=d_t_work(ind+j) + v(j)=incr(j)+d_t(j,inres) + enddo + else + do j=1,3 + v(j)=incr(j) enddo - ind=ind+3 endif - enddo + call vecpr(pr(1),v(1),vp) +!c write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3), +!c & " v",(v(j),j=1,3)," vp",(vp(j),j=1,3) +! if (mnum.gt.4) then +! mscab=0.0d0 +! else + mscab=msc(iti,mnum) +! endif + do j=1,3 + L(j)=L(j)+mscab*vp(j) + enddo +!c write (iout,*) "L",(l(j),j=1,3) + if (itype(i,1).ne.10 .and.itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.mnum.le.2) then + do j=1,3 + v(j)=incr(j)+d_t(j,inres) + enddo + call vecpr(dc(1,inres),d_t(1,inres),vp) + do j=1,3 + L(j)=L(j)+Isc(iti,mnum)*vp(j) + enddo + endif + do j=1,3 + incr(j)=incr(j)+d_t(j,i) + enddo + enddo return - end subroutine sd_verlet2_ciccotti -#endif -!----------------------------------------------------------------------------- -! moments.f -!----------------------------------------------------------------------------- + end subroutine angmom +!--------------------------------------------------------------- + subroutine vcm_vel(vcm) + double precision vcm(3),vv(3),summas,amas + integer i,j,mnum,mnum1 + do j=1,3 + vcm(j)=0.0d0 + vv(j)=d_t(j,0) + enddo + summas=0.0d0 + do i=nnt,nct + mnum=molnum(i) + if ((mnum.ge.5).or.(mnum.eq.3))& + mp(mnum)=msc(itype(i,mnum),mnum) + if (i.lt.nct) then + summas=summas+mp(mnum) + do j=1,3 + vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i)) + enddo + endif + if (mnum.ne.4) then + amas=msc(iabs(itype(i,mnum)),mnum) + else + amas=0.0d0 + endif +! amas=msc(iabs(itype(i))) + summas=summas+amas +! if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then + if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.(mnum.lt.4)) then + do j=1,3 + vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres)) + enddo + else + do j=1,3 + vcm(j)=vcm(j)+amas*vv(j) + enddo + endif + do j=1,3 + vv(j)=vv(j)+d_t(j,i) + enddo + enddo +!c write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas + do j=1,3 + vcm(j)=vcm(j)/summas + enddo + return + end subroutine vcm_vel +#else subroutine inertia_tensor ! Calculating the intertia tensor for the entire protein in order to @@ -3878,9 +4645,11 @@ M_PEP=0.0d0 do i=nnt,nct-1 mnum=molnum(i) - if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum) - if (itype(i,mnum).eq.ntyp1_molec(mnum)) cycle + if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum) + write(iout,*) "WTF",itype(i,mnum),i,mnum,mp(mnum) +! if (itype(i,mnum).eq.ntyp1_molec(mnum)) cycle M_PEP=M_PEP+mp(mnum) + do j=1,3 cm(j)=cm(j)+(c(j,i)+0.5d0*dc(j,i))*mp(mnum) enddo @@ -3891,7 +4660,7 @@ M_SC=0.0d0 do i=nnt,nct mnum=molnum(i) - if (mnum.eq.5) cycle + if (mnum.ge.5) cycle iti=iabs(itype(i,mnum)) M_SC=M_SC+msc(iabs(iti),mnum) inres=i+nres @@ -3902,10 +4671,10 @@ do j=1,3 cm(j)=cm(j)/(M_SC+M_PEP) enddo - +! write(iout,*) "Center of mass:",cm do i=nnt,nct-1 mnum=molnum(i) - if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum) + if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum) do j=1,3 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j) enddo @@ -3916,11 +4685,16 @@ Im(2,2)=Im(2,2)+mp(mnum)*(pr(3)*pr(3)+pr(1)*pr(1)) Im(3,3)=Im(3,3)+mp(mnum)*(pr(1)*pr(1)+pr(2)*pr(2)) enddo + +! write(iout,*) "The angular momentum before msc add" +! do i=1,3 +! write (iout,*) (Im(i,j),j=1,3) +! enddo do i=nnt,nct mnum=molnum(i) iti=iabs(itype(i,mnum)) - if (mnum.eq.5) cycle + if (mnum.ge.5) cycle inres=i+nres do j=1,3 pr(j)=c(j,inres)-cm(j) @@ -3932,6 +4706,10 @@ Im(2,2)=Im(2,2)+msc(iabs(iti),mnum)*(pr(3)*pr(3)+pr(1)*pr(1)) Im(3,3)=Im(3,3)+msc(iabs(iti),mnum)*(pr(1)*pr(1)+pr(2)*pr(2)) enddo +! write(iout,*) "The angular momentum before Ip add" +! do i=1,3 +! write (iout,*) (Im(i,j),j=1,3) +! enddo do i=nnt,nct-1 mnum=molnum(i) @@ -3948,13 +4726,17 @@ Im(3,3)=Im(3,3)+Ip(mnum)*(1-dc_norm(3,i)*dc_norm(3,i))* & vbld(i+1)*vbld(i+1)*0.25d0 enddo +! write(iout,*) "The angular momentum before Isc add" +! do i=1,3 +! write (iout,*) (Im(i,j),j=1,3) +! enddo do i=nnt,nct mnum=molnum(i) ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& - .and.(mnum.ne.5)) then + .and.(mnum.lt.4)) then iti=iabs(itype(i,mnum)) inres=i+nres Im(1,1)=Im(1,1)+Isc(iti,mnum)*(1-dc_norm(1,inres)* & @@ -3972,10 +4754,17 @@ endif enddo +! write(iout,*) "The angular momentum before agnom:" +! do i=1,3 +! write (iout,*) (Im(i,j),j=1,3) +! enddo + call angmom(cm,L) ! write(iout,*) "The angular momentum before adjustment:" ! write(iout,*) (L(j),j=1,3) - +! do i=1,3 +! write (iout,*) (Im(i,j),j=1,3) +! enddo Im(2,1)=Im(1,2) Im(3,1)=Im(1,3) Im(3,2)=Im(2,3) @@ -4033,13 +4822,15 @@ do i=nnt,nct-1 call vecpr(vrot(1),dc(1,i),vp) do j=1,3 +! print *,"HERE2",d_t(j,i),vp(j) d_t(j,i)=d_t(j,i)-vp(j) +! print *,"HERE2",d_t(j,i),vp(j) enddo enddo do i=nnt,nct mnum=molnum(i) if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& - .and.(mnum.ne.5)) then + .and.(mnum.lt.4)) then ! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then ! if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then inres=i+nres @@ -4083,7 +4874,7 @@ enddo do i=nnt,nct-1 mnum=molnum(i) - if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum) + if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum) do j=1,3 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j) enddo @@ -4096,12 +4887,14 @@ call vecpr(pr(1),v(1),vp) do j=1,3 L(j)=L(j)+mp(mnum)*vp(j) +! print *,"HERE3",J,i,L(j),mp(mnum),Ip(mnum),mnum enddo do j=1,3 pr(j)=0.5d0*dc(j,i) pp(j)=0.5d0*d_t(j,i) enddo call vecpr(pr(1),pp(1),vp) +! print *,vp,"vp" do j=1,3 L(j)=L(j)+Ip(mnum)*vp(j) enddo @@ -4113,7 +4906,7 @@ mnum=molnum(i) iti=iabs(itype(i,mnum)) inres=i+nres - if (mnum.eq.5) then + if (mnum.gt.4) then mscab=0.0d0 else mscab=msc(iti,mnum) @@ -4121,8 +4914,9 @@ do j=1,3 pr(j)=c(j,inres)-cm(j) enddo + !endif if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& - .and.(mnum.ne.5)) then + .and.(mnum.lt.4)) then do j=1,3 v(j)=incr(j)+d_t(j,inres) enddo @@ -4131,6 +4925,7 @@ v(j)=incr(j) enddo endif +! print *,i,pr,"pr",v call vecpr(pr(1),v(1),vp) ! write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),& ! " v",(v(j),j=1,3)," vp",(vp(j),j=1,3) @@ -4138,8 +4933,10 @@ L(j)=L(j)+mscab*vp(j) enddo ! write (iout,*) "L",(l(j),j=1,3) +! print *,"L",(l(j),j=1,3),i,vp(1) + if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& - .and.(mnum.ne.5)) then + .and.(mnum.lt.4)) then do j=1,3 v(j)=incr(j)+d_t(j,inres) enddo @@ -4179,21 +4976,22 @@ summas=0.0d0 do i=nnt,nct mnum=molnum(i) - if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum) + if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum) if (i.lt.nct) then summas=summas+mp(mnum) do j=1,3 vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i)) +! print *,i,j,vv(j),d_t(j,i) enddo endif - if (mnum.ne.5) then + if (mnum.ne.4) then amas=msc(iabs(itype(i,mnum)),mnum) else amas=0.0d0 endif summas=summas+amas if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& - .and.(mnum.ne.5)) then + .and.(mnum.lt.4)) then do j=1,3 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres)) enddo @@ -4212,6 +5010,7 @@ enddo return end subroutine vcm_vel +#endif !----------------------------------------------------------------------------- ! rattle.F !----------------------------------------------------------------------------- @@ -4264,7 +5063,7 @@ do i=nnt,nct mnum=molnum(i) if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& - .and.(mnum.ne.5)) nbond=nbond+1 + .and.(mnum.lt.4)) nbond=nbond+1 enddo ! Make a folded form of the Ginv-matrix ind=0 @@ -4285,7 +5084,7 @@ do k=nnt,nct mnum=molnum(k) if (itype(k,1).ne.10 .and. itype(k,mnum).ne.ntyp1_molec(mnum)& - .and.(mnum.ne.5)) then + .and.(mnum.lt.4)) then jj=jj+1 do l=1,3 ind1=ind1+1 @@ -4298,7 +5097,7 @@ do i=nnt,nct mnum=molnum(i) if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& - .and.(mnum.ne.5)) + .and.(mnum.lt.4)) ii=ii+1 do j=1,3 ind=ind+1 @@ -4396,7 +5195,7 @@ do i=nnt,nct mnum=molnum(i) if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& - .and.(mnum.ne.5)) then + .and.(mnum.lt.4)) then ind=ind+1 write (iout,'(2i5,3f10.5,5x,e15.5)') & @@ -4479,7 +5278,7 @@ do i=nnt,nct mnum=molnum(i) if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& - .and.(mnum.ne.5)) + .and.(mnum.lt.4)) ind=ind+1 xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') & @@ -4571,7 +5370,7 @@ do k=nnt,nct mnum=molnum(i) if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& - .and.(mnum.ne.5)) then + .and.(mnum.lt.4)) then ind=ind+1 do j=1,3 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k+nres) @@ -4602,7 +5401,7 @@ do i=nnt,nct mnum=molnum(i) if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& - .and.(mnum.ne.5)) then + .and.(mnum.lt.4)) then ind=ind+1 do j=1,nbond Cmat(ind,j)=0.0d0 @@ -4621,7 +5420,7 @@ do i=nnt,nct mnum=molnum(i) if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& - .and.(mnum.ne.5)) then + .and.(mnum.lt.4)) then ind=ind+1 x(ind)=scalar(d_t(1,i+nres),dC(1,i+nres)) endif @@ -4637,7 +5436,7 @@ do i=nnt,nct mnum=molnum(i) if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& - .and.(mnum.ne.5)) then + .and.(mnum.lt.4)) then ind=ind+1 write (iout,'(2i5,3f10.5,5x,e15.5)') & i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind) @@ -4674,7 +5473,7 @@ do i=nnt,nct mnum=molnum(i) if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& - .and.(mnum.ne.5)) then + .and.(mnum.lt.4)) then ind=ind+1 do j=1,3 xx=0.0d0 @@ -4697,7 +5496,7 @@ do i=nnt,nct mnum=molnum(i) if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& - .and.(mnum.ne.5)) + .and.(mnum.lt.4)) ind=ind+1 write (iout,'(2i5,3f10.5,5x,2e15.5)') & i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind),& @@ -5029,15 +5828,113 @@ ! include 'COMMON.IOUNITS' !el real(kind=8),dimension(6*nres) :: gamvec !(MAXRES6) maxres6=6*maxres !el common /syfek/ gamvec - real(kind=8) :: vv(3),vvtot(3,nres),v_work(6*nres) !,& +#ifdef FIVEDIAG + integer iposc,ichain,n,innt,inct + double precision rs(nres*2) + real(kind=8) ::v_work(3,6*nres),vvec(2*nres) +#else + real(kind=8) :: v_work(6*nres) +#endif + + real(kind=8) :: vv(3),vvtot(3,nres)!,v_work(6*nres) !,& !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. +#ifdef FIVEDIAG +!c Here accelerations due to friction forces are computed right after forces. + d_t_work(:6*nres)=0.0d0 + do j=1,3 + v_work(j,1)=d_t(j,0) + v_work(j,nnt)=d_t(j,0) + enddo + do i=nnt+1,nct + do j=1,3 + v_work(j,i)=v_work(j,i-1)+d_t(j,i-1) + enddo + enddo + do i=nnt,nct + mnum=molnum(i) + if (iabs(itype(i,1)).ne.10 .and. iabs(itype(i,mnum)).ne.ntyp1_molec(mnum).and.mnum.lt.3) then + do j=1,3 + v_work(j,i+nres)=v_work(j,i)+d_t(j,i+nres) + enddo + endif + enddo +#ifdef DEBUG + write (iout,*) "v_work" + do i=1,2*nres + write (iout,'(i5,3f10.5)') i,(v_work(j,i),j=1,3) + enddo +#endif + do j=1,3 + ind=0 + do ichain=1,nchain + n=dimen_chain(ichain) + iposc=iposd_chain(ichain) +!c write (iout,*) "friction_force j",j," ichain",ichain, +!c & " n",n," iposc",iposc,iposc+n-1 + innt=chain_border(1,ichain) + inct=chain_border(2,ichain) +!c diagnostics +!c innt=chain_border(1,1) +!c inct=chain_border(2,1) + do i=innt,inct + vvec(ind+1)=v_work(j,i) + ind=ind+1 +! if (iabs(itype(i)).ne.10) then + if (iabs(itype(i,1)).ne.10 .and. iabs(itype(i,mnum)).ne.ntyp1_molec(mnum).and.mnum.lt.3) then + vvec(ind+1)=v_work(j,i+nres) + ind=ind+1 + endif + enddo +#ifdef DEBUG + write (iout,*) "vvec ind",ind," n",n + write (iout,'(f10.5)') (vvec(i),i=iposc,ind) +#endif +!c write (iout,*) "chain",i," ind",ind," n",n +#ifdef TIMING +#ifdef MPI + time01=MPI_Wtime() +#else + time01=tcpu() +#endif +#endif +! if (large) print *,"before fivediagmult" + call fivediagmult(n,DMfric(iposc),DU1fric(iposc),& + DU2fric(iposc),vvec(iposc),rs) +! if (large) print *,"after fivediagmult" +#ifdef TIMING +#ifdef MPI + time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01 +#else + time_fricmatmult=time_fricmatmult+tcpu()-time01 +#endif +#endif +#ifdef DEBUG + write (iout,*) "rs" + write (iout,'(f10.5)') (rs(i),i=1,n) +#endif + do i=iposc,iposc+n-1 +! write (iout,*) "ichain",ichain," i",i," j",j,& +! "index",3*(i-1)+j,"rs",rs(i-iposc+1) + fric_work(3*(i-1)+j)=-rs(i-iposc+1) + enddo + enddo + enddo +#ifdef DEBUG + write (iout,*) "Vector fric_work dimen3",dimen3 + write (iout,'(3f10.5)') (fric_work(j),j=1,dimen3) +#endif +#else if(.not.allocated(gamvec)) allocate(gamvec(nres6)) !(MAXRES6) if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres do i=0,nres2 @@ -5059,7 +5956,7 @@ do i=nnt,nct mnum=molnum(i) if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))& - .and.(mnum.ne.5)) then + .and.(mnum.lt.4)) then do j=1,3 d_t_work(ind+j)=d_t(j,i+nres) enddo @@ -5090,7 +5987,7 @@ do i=nnt,nct mnum=molnum(i) if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))& - .and.(mnum.ne.5)) then + .and.(mnum.lt.4)) then ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then do j=1,3 friction(j,i+nres)=fric_work(ind+j) @@ -5173,6 +6070,7 @@ enddo enddo endif +#endif return end subroutine friction_force !----------------------------------------------------------------------------- @@ -5208,7 +6106,7 @@ !#endif ! include 'COMMON.IOUNITS' integer :: IERROR - integer :: i,j,ind,ind1,m + integer :: i,j,ind,ind1,m,ichain,innt,inct logical :: lprn = .false. real(kind=8) :: dtdi !el ,gamvec(2*nres) !el real(kind=8),dimension(2*nres,2*nres) :: ginvfric,fcopy @@ -5222,20 +6120,155 @@ nres2=2*nres nres6=6*nres #ifdef MPI +#ifndef FIVEDIAG if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2)) if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres if (fg_rank.ne.king) goto 10 #endif +#endif ! nres2=2*nres ! nres6=6*nres if(.not.allocated(gamvec)) allocate(gamvec(nres2)) !(MAXRES2) +#ifndef FIVEDIAG if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres !el allocate(fcopy(nres2,nres2)) !maxres2=2*maxres if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !maxres2=2*maxres if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2)) +#endif +#ifdef FIVEDIAG + if (.not.allocated(DMfric)) allocate(DMfric(nres2)) + if (.not.allocated(DU1fric)) allocate(DU1fric(nres2)) + if (.not.allocated(DU2fric)) allocate(DU2fric(nres2)) +! Load the friction coefficients corresponding to peptide groups + ind1=0 + do i=nnt,nct-1 + mnum=molnum(i) + ind1=ind1+1 + gamvec(ind1)=gamp(mnum) + enddo +!HERE TEST + if (molnum(nct).eq.5) then + mnum=molnum(i) + ind1=ind1+1 + gamvec(ind1)=gamp(mnum) + endif +! Load the friction coefficients corresponding to side chains + m=nct-nnt + ind=0 + do j=1,2 + gamsc(ntyp1_molec(j),j)=1.0d0 + enddo + do i=nnt,nct + mnum=molnum(i) + ind=ind+1 + ii = ind+m + iti=itype(i,mnum) + gamvec(ii)=gamsc(iabs(iti),mnum) + enddo + if (surfarea) call sdarea(gamvec) + DMfric=0.0d0 + DU1fric=0.0d0 + DU2fric=0.0d0 + ind=1 + do ichain=1,nchain + innt=chain_border(1,ichain) + inct=chain_border(2,ichain) +!c write (iout,*) "ichain",ichain," innt",innt," inct",inct +!c DMfric part + mnum=molnum(innt) + DMfric(ind)=gamvec(innt-nnt+1)/4 + if (iabs(itype(innt,1)).eq.10.or.mnum.gt.2) then + DMfric(ind)=DMfric(ind)+gamvec(m+innt-nnt+1) + ind=ind+1 + else + DMfric(ind+1)=gamvec(m+innt-nnt+1) + ind=ind+2 + endif +!c write (iout,*) "DMfric init ind",ind +!c DMfric + do i=innt+1,inct-1 + mnum=molnum(i) + DMfric(ind)=gamvec(i-nnt+1)/2 + if (iabs(itype(i,1)).eq.10.or.mnum.gt.2) then + DMfric(ind)=DMfric(ind)+gamvec(m+i-nnt+1) + ind=ind+1 + else + DMfric(ind+1)=gamvec(m+i-nnt+1) + ind=ind+2 + endif + enddo +!c write (iout,*) "DMfric endloop ind",ind + if (inct.gt.innt) then + DMfric(ind)=gamvec(inct-1-nnt+1)/4 + mnum=molnum(inct) + if (iabs(itype(inct,1)).eq.10.or.mnum.gt.2) then + DMfric(ind)=DMfric(ind)+gamvec(inct+m-nnt+1) + ind=ind+1 + else + DMfric(ind+1)=gamvec(inct+m-nnt+1) + ind=ind+2 + endif + endif +!c write (iout,*) "DMfric end ind",ind + enddo +!c DU1fric part + do ichain=1,nchain + mnum=molnum(i) + + ind=iposd_chain(ichain) + innt=chain_border(1,ichain) + inct=chain_border(2,ichain) + do i=innt,inct + if (iabs(itype(i,1)).ne.10.and.mnum.le.2) then + ind=ind+2 + else + DU1fric(ind)=gamvec(i-nnt+1)/4 + ind=ind+1 + endif + enddo + enddo +!c DU2fric part + do ichain=1,nchain + mnum=molnum(i) + ind=iposd_chain(ichain) + innt=chain_border(1,ichain) + inct=chain_border(2,ichain) + do i=innt,inct-1 + if (iabs(itype(i,1)).ne.10.and.mnum.le.2) then + DU2fric(ind)=gamvec(i-nnt+1)/4 + DU2fric(ind+1)=0.0d0 + ind=ind+2 + else + DU2fric(ind)=0.0d0 + ind=ind+1 + endif + enddo + enddo + if (lprn) then + write(iout,*)"The upper part of the five-diagonal friction matrix" + do ichain=1,nchain + write (iout,'(a,i5)') 'Chain',ichain + innt=iposd_chain(ichain) + inct=iposd_chain(ichain)+dimen_chain(ichain)-1 + do i=innt,inct + if (i.lt.inct-1) then + write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i),DU1fric(i),& + DU2fric(i) + else if (i.eq.inct-1) then + write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i),DU1fric(i) + else + write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i) + endif + enddo + enddo + endif + 10 continue +#else + + ! Zeroing out fricmat do i=1,dimen do j=1,dimen @@ -5249,6 +6282,12 @@ ind1=ind1+1 gamvec(ind1)=gamp(mnum) enddo +!HERE TEST + if (molnum(nct).eq.5) then + mnum=molnum(i) + ind1=ind1+1 + gamvec(ind1)=gamp(mnum) + endif ! Load the friction coefficients corresponding to side chains m=nct-nnt ind=0 @@ -5421,6 +6460,7 @@ ! call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy) endif #endif +#endif return end subroutine setup_fricmat !----------------------------------------------------------------------------- @@ -5456,6 +6496,9 @@ real(kind=8) :: time00 logical :: lprn = .false. integer :: i,j,ind,mnum +#ifdef FIVEDIAG + integer ichain,innt,inct,iposc +#endif do i=0,2*nres do j=1,3 @@ -5491,6 +6534,72 @@ #else time_fsample=time_fsample+tcpu()-time00 #endif +#ifdef FIVEDIAG + ind=0 + do ichain=1,nchain + innt=chain_border(1,ichain) + inct=chain_border(2,ichain) + iposc=iposd_chain(ichain) +!c for debugging only +!c innt=chain_border(1,1) +!c inct=chain_border(2,1) +!c iposc=iposd_chain(1) +!c write (iout,*)"stochastic_force ichain=",ichain," innt",innt, +!c & " inct",inct," iposc",iposc + do j=1,3 + stochforcvec(ind+j)=0.5d0*force(j,innt) + enddo + if (iabs(itype(innt,molnum(innt))).eq.10.or.molnum(innt).ge.3) then + do j=1,3 + stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,innt+nres) + enddo + ind=ind+3 + else + ind=ind+3 + do j=1,3 + stochforcvec(ind+j)=force(j,innt+nres) + enddo + ind=ind+3 + endif + do i=innt+1,inct-1 + do j=1,3 + stochforcvec(ind+j)=0.5d0*(force(j,i)+force(j,i-1)) + enddo + if (iabs(itype(i,1).eq.10).or.molnum(i).ge.3) then + do j=1,3 + stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,i+nres) + enddo + ind=ind+3 + else + ind=ind+3 + do j=1,3 + stochforcvec(ind+j)=force(j,i+nres) + enddo + ind=ind+3 + endif + enddo + do j=1,3 + stochforcvec(ind+j)=0.5d0*force(j,inct-1) + enddo + if (iabs(itype(inct,molnum(inct))).eq.10.or.molnum(inct).ge.3) then + do j=1,3 + stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,inct+nres) + enddo + ind=ind+3 + else + ind=ind+3 + do j=1,3 + stochforcvec(ind+j)=force(j,inct+nres) + enddo + ind=ind+3 + endif +!c write (iout,*) "chain",ichain," ind",ind + enddo +#ifdef DEBUG + write (iout,*) "stochforcvec" + write (iout,'(3f10.5)') (stochforcvec(j),j=1,ind) +#endif +#else ! Compute the stochastic forces acting on virtual-bond vectors. do j=1,3 ff(j)=0.0d0 @@ -5517,7 +6626,7 @@ do i=nnt,nct mnum=molnum(i) if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))& - .and.(mnum.ne.5)) then + .and.(mnum.lt.4)) then ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then do j=1,3 stochforc(j,i+nres)=force(j,i+nres) @@ -5538,7 +6647,7 @@ do i=nnt,nct mnum=molnum(i) if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))& - .and.(mnum.ne.5)) then + .and.(mnum.lt.4)) then ! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then do j=1,3 stochforcvec(ind+j)=stochforc(j,i+nres) @@ -5598,7 +6707,7 @@ enddo endif - +#endif return end subroutine stochastic_force !----------------------------------------------------------------------------- @@ -6237,13 +7346,16 @@ ! commom.langevin.lang0 ! common /langforc/ allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2) +#ifndef FIVEDIAG if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2)) allocate(fricvec(nres2,nres2)) !(MAXRES2,MAXRES2) +#endif allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6) allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch) #endif - +#ifndef FIVEDIAG if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) +#endif !---------------------- ! commom.hairpin in CSA module !---------------------- @@ -6262,6 +7374,9 @@ #ifdef FIVEDIAG allocate(DM(nres2),DU1(nres2),DU2(nres2)) allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2)) +!#ifdef DEBUG + allocate(Gvec(nres2,nres2)) +!#endif #else write (iout,*) "Before A Allocation",nfgtasks-1 call flush(iout)