X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2FMD.f90;h=e34b3840191dabba240fbe66f40af5d196af8eb6;hb=10689ab7d813dfbdbb0c6e631d90234d78ea306a;hp=edbd8102d2425edb107b2633042c36479fff787d;hpb=ae71f1947b5c88582ff30337921f4bee3d5d3038;p=unres4.git

diff --git a/source/unres/MD.f90 b/source/unres/MD.f90
index edbd810..e34b384 100644
--- a/source/unres/MD.f90
+++ b/source/unres/MD.f90
@@ -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
 
@@ -978,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)
@@ -1013,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)
@@ -2493,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)
@@ -2693,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
@@ -2952,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
@@ -2977,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
@@ -3002,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
@@ -3065,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
@@ -3855,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)
@@ -3865,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
          	   		  			      
@@ -3932,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
@@ -4013,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