X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2FMD.f90;h=e34b3840191dabba240fbe66f40af5d196af8eb6;hb=10689ab7d813dfbdbb0c6e631d90234d78ea306a;hp=1ddffac5517567bd1fb35588afe8177ab364f21e;hpb=1b595d5c963c2f521e7b7e6c426dcc358feaf430;p=unres4.git diff --git a/source/unres/MD.f90 b/source/unres/MD.f90 index 1ddffac..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 @@ -832,7 +841,7 @@ #else tt0 = tcpu() #endif - print *,"just befor setup matix" + print *,"just befor setup matix",nres ! Determine the inverse of the inertia matrix. call setup_MD_matrices ! Initialize MD @@ -883,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 @@ -976,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) @@ -1011,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) @@ -2491,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) @@ -2691,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 @@ -2950,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 @@ -2975,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 @@ -3000,7 +3018,7 @@ ! 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)) + .or.(mnum.eq.5)) then do j=1,3 d_t(j,i)=d_t(j,i+1)-d_t(j,i) enddo @@ -3063,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 @@ -3853,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) @@ -3863,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 @@ -3930,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 @@ -4011,8 +4029,8 @@ 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)+mscab*vp(j) enddo @@ -5091,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 @@ -6118,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 !----------------------