X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2FMD.f90;h=32e48f80e78195343b8869db9e960f8f15baed3b;hb=a30537f5834fb1220b66036bc509bb3ab46dac7a;hp=848702bbc90cdf80b147b465faa46b5372875420;hpb=2dfae3f5bcdd236ada1c317c08de18027bda3b7d;p=unres4.git diff --git a/source/unres/MD.f90 b/source/unres/MD.f90 index 848702b..32e48f8 100644 --- a/source/unres/MD.f90 +++ b/source/unres/MD.f90 @@ -141,7 +141,7 @@ real(kind=8),dimension(:,:),allocatable :: ginvfric !(2*nres,2*nres) !maxres2=2*maxres !----------------------------------------------------------------------------- ! common /przechowalnia/ subroutine: setup_fricmat -!el real(kind=8),dimension(:,:),allocatable :: fcopy !(2*nres,2*nres) + real(kind=8),dimension(:,:),allocatable :: fcopy !(2*nres,2*nres) !----------------------------------------------------------------------------- ! ! @@ -652,7 +652,7 @@ ! include 'COMMON.INTERACT' ! include 'COMMON.MD' ! include 'COMMON.IOUNITS' - real(kind=8) :: KE_total + real(kind=8) :: KE_total,mscab integer :: i,j,k,iti,mnum real(kind=8) :: KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),& @@ -673,6 +673,7 @@ enddo do i=nnt,nct-1 mnum=molnum(i) + if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum) ! write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3) do j=1,3 v(j)=incr(j)+0.5d0*d_t(j,i) @@ -693,9 +694,14 @@ do i=nnt,nct mnum=molnum(i) iti=iabs(itype(i,mnum)) + if (mnum.eq.5) then + mscab=0.0d0 + else + mscab=msc(iti,mnum) + endif ! 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 + if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)& + .or.(mnum.eq.5)) then do j=1,3 v(j)=incr(j) enddo @@ -706,7 +712,7 @@ endif ! write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3) ! 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)) + 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 incr(j)=incr(j)+d_t(j,i) @@ -804,8 +810,10 @@ real(kind=8) :: tt0,scalfac integer :: nres2 nres2=2*nres + print *, "ENTER MD" ! #ifdef MPI + print *,"MY tmpdir",tmpdir,ilen(tmpdir) if (ilen(tmpdir).gt.0) & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_" & //liczba(:ilen(liczba))//'.rst') @@ -824,10 +832,14 @@ #else tt0 = tcpu() #endif + print *,"just befor setup matix",nres ! Determine the inverse of the inertia matrix. call setup_MD_matrices ! Initialize MD + print *,"AFTER SETUP MATRICES" call init_MD + print *,"AFTER INIT MD" + #ifdef MPI t_MDsetup = MPI_Wtime()-tt0 #else @@ -871,7 +883,9 @@ stop #endif else if (lang.eq.1 .or. lang.eq.4) then + print *,"before setup_fricmat" call setup_fricmat + print *,"after setup_fricmat" endif #ifdef MPI t_langsetup=MPI_Wtime()-tt0 @@ -964,7 +978,10 @@ #endif endif if (ntwe.ne.0) then - if (mod(itime,ntwe).eq.0) call statout(itime) + if (mod(itime,ntwe).eq.0) then + call statout(itime) + call returnbox + endif #ifdef VOUT do j=1,3 v_work(j)=d_t(j,0) @@ -977,7 +994,8 @@ enddo enddo do i=nnt,nct - if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then + mnum=molnum(i) + if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.mnum.ne.5) then do j=1,3 ind=ind+1 v_work(ind)=d_t(j,i+nres) @@ -998,6 +1016,7 @@ #endif endif if (mod(itime,ntwx).eq.0) then + call returnbox write (tytul,'("time",f8.2)') totT if(mdpdb) then call hairpin(.true.,nharp,iharp) @@ -2209,7 +2228,8 @@ enddo endif ! Side chains - if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then + if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.& + molnum(i).ne.5) then do j=1,3 epdriftij= & dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i)) @@ -2911,10 +2931,11 @@ 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).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& - .and.(mnum.ne.5)) 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 @@ -2984,8 +3005,8 @@ do i=nnt,nct ! if (itype(i,1).eq.10) then 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).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)& + .or.(mnum.eq.5)) do j=1,3 d_t(j,i)=d_t(j,i+1)-d_t(j,i) enddo @@ -3742,6 +3763,7 @@ 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 M_PEP=M_PEP+mp(mnum) do j=1,3 @@ -3754,6 +3776,7 @@ M_SC=0.0d0 do i=nnt,nct mnum=molnum(i) + if (mnum.eq.5) cycle iti=iabs(itype(i,mnum)) M_SC=M_SC+msc(iabs(iti),mnum) inres=i+nres @@ -3767,6 +3790,7 @@ do i=nnt,nct-1 mnum=molnum(i) + if (mnum.eq.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 @@ -3781,6 +3805,7 @@ do i=nnt,nct mnum=molnum(i) iti=iabs(itype(i,mnum)) + if (mnum.eq.5) cycle inres=i+nres do j=1,3 pr(j)=c(j,inres)-cm(j) @@ -3931,7 +3956,7 @@ ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' - + real(kind=8) :: mscab real(kind=8),dimension(3) :: L,cm,pr,vp,vrot,incr,v,pp integer :: iti,inres,i,j,mnum ! Calculate the angular momentum @@ -3943,6 +3968,7 @@ enddo do i=nnt,nct-1 mnum=molnum(i) + if (mnum.eq.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 @@ -3972,6 +3998,11 @@ mnum=molnum(i) iti=iabs(itype(i,mnum)) inres=i+nres + if (mnum.eq.5) then + mscab=0.0d0 + else + mscab=msc(iti,mnum) + endif do j=1,3 pr(j)=c(j,inres)-cm(j) enddo @@ -3989,7 +4020,7 @@ ! write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3), ! & " v",(v(j),j=1,3)," vp",(vp(j),j=1,3) do j=1,3 - L(j)=L(j)+msc(iabs(iti),mnum)*vp(j) + L(j)=L(j)+mscab*vp(j) enddo ! write (iout,*) "L",(l(j),j=1,3) if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& @@ -4033,13 +4064,18 @@ summas=0.0d0 do i=nnt,nct mnum=molnum(i) + if (mnum.eq.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)) enddo endif + if (mnum.ne.5) 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 @@ -5061,26 +5097,30 @@ logical :: lprn = .false. real(kind=8) :: dtdi !el ,gamvec(2*nres) !el real(kind=8),dimension(2*nres,2*nres) :: ginvfric,fcopy - real(kind=8),dimension(2*nres,2*nres) :: fcopy +! real(kind=8),allocatable,dimension(:,:) :: fcopy !el real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2)) !el common /syfek/ gamvec real(kind=8) :: work(8*2*nres) integer :: iwork(2*nres) !el common /przechowalnia/ ginvfric,Ghalf,fcopy integer :: ii,iti,k,l,nzero,nres2,nres6,ierr,mnum + nres2=2*nres + nres6=6*nres #ifdef MPI + 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 - nres2=2*nres - nres6=6*nres +! nres2=2*nres +! nres6=6*nres if(.not.allocated(gamvec)) allocate(gamvec(nres2)) !(MAXRES2) if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres -!el if(.not.allocated(fcopy)) allocate(fcopy(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 -!el if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2)) + if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2)) ! Zeroing out fricmat do i=1,dimen do j=1,dimen @@ -5097,7 +5137,7 @@ ! Load the friction coefficients corresponding to side chains m=nct-nnt ind=0 - do j=1,5 + do j=1,2 gamsc(ntyp1_molec(j),j)=1.0d0 enddo do i=nnt,nct @@ -6088,7 +6128,7 @@ allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch) #endif -!el if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) + if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !---------------------- ! commom.hairpin in CSA module !----------------------