X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2FMD.f90;h=e34b3840191dabba240fbe66f40af5d196af8eb6;hb=10689ab7d813dfbdbb0c6e631d90234d78ea306a;hp=848702bbc90cdf80b147b465faa46b5372875420;hpb=2dfae3f5bcdd236ada1c317c08de18027bda3b7d;p=unres4.git diff --git a/source/unres/MD.f90 b/source/unres/MD.f90 index 848702b..e34b384 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) !----------------------------------------------------------------------------- ! ! @@ -187,6 +187,9 @@ real(kind=8),dimension(6*nres,2*nres) :: Bmat,GBmat,Tmat !(MAXRES6,MAXRES2) (maxres2=2*maxres,maxres6=6*maxres) real(kind=8),dimension(2*nres,2*nres) :: Cmat_,Cinv !(maxres2,maxres2) maxres2=2*maxres real(kind=8),dimension(6*nres,6*nres) :: Pmat !(maxres6,maxres6) maxres6=6*maxres +! real(kind=8),dimension(:,:),allocatable :: Bmat,GBmat,Tmat !(MAXRES6,MAXRES2) (maxres2=2*maxres,maxres6=6*maxres) +! real(kind=8),dimension(:,:),allocatable :: Cmat_,Cinv !(maxres2,maxres2) maxres2=2*maxres +! real(kind=8),dimension(:,:),allocatable :: Pmat !(maxres6,maxres6) maxres6=6*maxres real(kind=8),dimension(6*nres) :: Td !(maxres6) maxres6=6*maxres real(kind=8),dimension(2*nres) :: ppvec !(maxres2) maxres2=2*maxres !el common /stochcalc/ stochforcvec @@ -202,6 +205,12 @@ logical :: osob nres2=2*nres nres6=6*nres +! if (.not.allocated(Bmat)) allocate(Bmat(nres6,nres2)) +! if (.not.allocated(GBmat)) allocate (GBmat(nres6,nres2)) +! if (.not.allocated(Tmat)) allocate (Tmat(nres6,nres2)) +! if (.not.allocated(Cmat_)) allocate(Cmat_(nres2,nres2)) +! if (.not.allocated(Cinv)) allocate (Cinv(nres2,nres2)) +! if (.not.allocated(Pmat)) allocate(Pmat(6*nres,6*nres)) if (.not.allocated(stochforcvec)) allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres @@ -652,7 +661,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 +682,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 +703,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 +721,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 +819,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 +841,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 +892,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 +987,11 @@ #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 +! call check_ecartint + endif #ifdef VOUT do j=1,3 v_work(j)=d_t(j,0) @@ -977,7 +1004,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 +1026,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 +2238,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)) @@ -2477,9 +2507,9 @@ write (iout,*) "vcm right after adjustment:" write (iout,*) (vcm(j),j=1,3) endif - if ((.not.rest).and.(indpdb.eq.0)) then - call chainbuild - if(iranconf.ne.0) then + if (.not.rest) then +! call chainbuild + if(iranconf.ne.0 .or.indpdb.gt.0.and..not.unres_pdb .or.preminim) then if (overlapsc) then print *, 'Calling OVERLAP_SC' call overlap_sc(fail) @@ -2677,12 +2707,14 @@ ! include 'COMMON.NAMES' ! include 'COMMON.TIME1' real(kind=8) :: xv,sigv,lowb,highb ,Ek1 +!#define DEBUG #ifdef FIVEDIAG real(kind=8) ,allocatable, dimension(:) :: DDU1,DDU2,DL2,DL1,xsolv,DML,rs real(kind=8) :: sumx #ifdef DEBUG real(kind=8) ,allocatable, dimension(:) :: rsold real (kind=8),allocatable,dimension(:,:) :: matold + integer :: iti #endif #endif integer :: i,j,ii,k,ind,mark,imark,mnum @@ -2911,10 +2943,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 @@ -2935,7 +2968,7 @@ 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*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 enddo write (iout,*) "i",i," ii",ii @@ -2960,9 +2993,9 @@ d_t(k,j)=d_t_work(ind) ind=ind+1 enddo - mnum=molnum(i) - if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& - .and.(mnum.ne.5)) + 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) ind=ind+1 @@ -2984,8 +3017,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)) then do j=1,3 d_t(j,i)=d_t(j,i+1)-d_t(j,i) enddo @@ -3048,8 +3081,8 @@ ! write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",& ! 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1 ! call flush(iout) +! write(iout,*) "end init MD" return -#undef DEBUG end subroutine random_vel !----------------------------------------------------------------------------- #ifndef LANG0 @@ -3742,6 +3775,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 +3788,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 +3802,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 +3817,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) @@ -3834,7 +3871,7 @@ call angmom(cm,L) ! write(iout,*) "The angular momentum before adjustment:" -! write(iout,*) (L(j),j=1,3) +! write(iout,*) (L(j),j=1,3) Im(2,1)=Im(1,2) Im(3,1)=Im(1,3) @@ -3844,7 +3881,7 @@ do i=1,3 do j=1,3 Imcp(i,j)=Im(i,j) - Id(i,j)=0.0d0 + Id(i,j)=0.0d0 enddo enddo @@ -3911,7 +3948,7 @@ enddo call angmom(cm,L) ! write(iout,*) "The angular momentum after adjustment:" -! write(iout,*) (L(j),j=1,3) +! write(iout,*) (L(j),j=1,3) return end subroutine inertia_tensor @@ -3931,7 +3968,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 +3980,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 +4010,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 @@ -3986,10 +4029,10 @@ enddo endif 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) +! 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 +4076,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 +5109,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 +5149,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 +6140,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 !----------------------