X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2FMD.F90;h=bf993ebfed1dd4737812463e5fae8ca8fc5a9ca6;hb=refs%2Fheads%2FUCGM;hp=87f739439c34f9cbe7cab8d2783e3b1dbc3a70da;hpb=91b9f78d94b96277537615722323ebe03cc0a014;p=unres4.git

diff --git a/source/unres/MD.F90 b/source/unres/MD.F90
index 87f7394..bf993eb 100644
--- a/source/unres/MD.F90
+++ b/source/unres/MD.F90
@@ -141,7 +141,9 @@
       real(kind=8),dimension(:,:),allocatable :: ginvfric !(2*nres,2*nres) !maxres2=2*maxres
 !-----------------------------------------------------------------------------
 !      common /przechowalnia/ subroutine: setup_fricmat
+#ifndef LBFGS
       real(kind=8),dimension(:,:),allocatable :: fcopy !(2*nres,2*nres)
+#endif
 !-----------------------------------------------------------------------------
 !
 !
@@ -697,6 +699,11 @@
        mnum=molnum(i)
         iti=iabs(itype(i,mnum))
         if (mnum.eq.5) iti=itype(i,mnum)
+!        if (itype(i,mnum).eq.ntyp1_molec(mnum)) then
+!         do j=1,3
+!          incr(j)=d_t(j,i)
+!         enddo
+!        endif
         if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
            .or.mnum.ge.3) then
           do j=1,3
@@ -724,7 +731,7 @@
          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)
+         if (mnum.eq.5) Ip(mnum)=0.0
 !        write (iout,*) "i",i
 !        write (iout,*) "i",i," mag1",mag1," mag2",mag2
          do j=1,3
@@ -920,8 +927,8 @@
 
       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)
+      if (mnum.eq.5) mp(mnum)=0.0d0
+!      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))
@@ -935,7 +942,11 @@
 !c to the velocities of the first Calpha.
       do i=innt,inct
         mnum=molnum(i)
+        if (mnum.eq.5) then
+        iti=itype(i,mnum)
+        else
         iti=iabs(itype(i,mnum))
+        endif
         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
@@ -959,7 +970,7 @@
          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)
+         if (mnum.eq.5) Ip(mnum)=0.0d0
 !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)&
@@ -970,7 +981,13 @@
 !c  The rotational part of the side chain virtual bond
        do i=innt,inct
          mnum=molnum(i)
-         iti=iabs(itype(i,mnum))
+!         iti=iabs(itype(i,mnum))
+        if (mnum.eq.5) then
+        iti=itype(i,mnum)
+        else
+        iti=iabs(itype(i,mnum))
+        endif
+
 !         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
@@ -1033,7 +1050,7 @@
 !      include 'COMMON.NAMES'
 !      include 'COMMON.TIME1'
 !      include 'COMMON.HAIRPIN'
-      real(kind=8),dimension(3) :: L,vcm
+      real(kind=8),dimension(3) :: L,vcm,boxx
 #ifdef VOUT
       real(kind=8),dimension(6*nres) :: v_work,v_transf	!(maxres6) maxres6=6*maxres
 #endif
@@ -1043,11 +1060,16 @@
 !el      common /gucio/ cm
       integer :: i,j,nharp
       integer,dimension(4,nres) :: iharp	!(4,nres/3)(4,maxres/3)
+  
 !      logical :: ovrtim
       real(kind=8) :: tt0,scalfac
       integer :: nres2,itime
       nres2=2*nres
       print *, "ENTER MD"
+      boxx(1)=boxxsize
+      boxx(2)=boxysize
+      boxx(3)=boxzsize
+
 !
 #ifdef MPI
       print *,"MY tmpdir",tmpdir,ilen(tmpdir)
@@ -1172,7 +1194,18 @@
         endif
         if (reset_vel .and. tbf .and. lang.eq.0 &
             .and. mod(itime,count_reset_vel).eq.0) then
+          !WARP WATER
+          do i=1,nres
+           if (molnum(i).eq.5) then
+             call to_box(c(1,i),c(2,i),c(3,i))
+             do j=1,3 
+             if (c(j,i).le.0) c(j,i)=c(j,i)+boxx(j)
+             enddo
+           endif
+!           write(iout,*) "COORD",c(1,i),c(2,i),c(3,i)
+          enddo
           call random_vel
+          
           write(iout,'(a,f20.2)') &
            "Velocities reset to random values, time",totT	
           do i=0,2*nres
@@ -1266,8 +1299,11 @@
 
           write (tytul,'("time",f8.2)') totT
           if(mdpdb) then
+             write(iout,*) "before hairpin"
              call hairpin(.true.,nharp,iharp)
+             write(iout,*) "before secondary"
              call secondary2(.true.)
+             write(iout,*) "before pdbout"
              call pdbout(potE,tytul,ipdb)
 !             call enerprint(potEcomp)
           else 
@@ -2254,6 +2290,10 @@
 !      call ginv_mult(fric_work, d_af_work)
 !      call ginv_mult(stochforcvec, d_as_work)
 #ifdef FIVEDIAG
+       write(iout,*) "forces before fivediaginv"
+      do i=1,dimen*3
+       write(iout,*) "fricwork",i,fric_work(i)
+      enddo
       call fivediaginv_mult(dimen,fric_work, d_af_work)
       call fivediaginv_mult(dimen,stochforcvec, d_as_work)
       if (large) then
@@ -2364,9 +2404,11 @@
 !      include 'COMMON.INTERACT'
 !      include 'COMMON.IOUNITS'
 !      include 'COMMON.NAMES'
-      real(kind=8),dimension(6*nres) :: stochforcvec,d_as_work1	!(MAXRES6) maxres6=6*maxres
+      real(kind=8),dimension(:),allocatable :: stochforcvec,d_as_work1	!(MAXRES6) maxres6=6*maxres
       real(kind=8) :: cos60 = 0.5d0, sin60 = 0.86602540378443864676d0
       integer :: i,j,ind,inres,mnum
+      if (.not.allocated(stochforcvec)) allocate(stochforcvec(6*nres))
+      if (.not.allocated(d_as_work1)) allocate(d_as_work1(6*nres))
 ! Revised 3/31/05 AL: correlation between random contributions to 
 ! position and velocity increments included.
 ! The correlation coefficients are calculated at low-friction limit.
@@ -2649,7 +2691,10 @@
       character(len=50) :: tytul
       logical :: file_exist
 !el      common /gucio/ cm
-      integer :: i,j,ipos,iq,iw,nft_sc,iretcode,nfun,ierr,mnum,itime
+      integer :: i,j,ipos,iq,iw,nft_sc,iretcode,ierr,mnum,itime
+#ifndef LBFGS 
+      integer :: nfun
+#endif
       real(kind=8) :: etot,tt0
       logical :: fail
 
@@ -2799,8 +2844,8 @@
       call inertia_tensor  
 !  Getting the potential energy and forces and velocities and accelerations
       call vcm_vel(vcm)
-!      write (iout,*) "velocity of the center of the mass:"
-!      write (iout,*) (vcm(j),j=1,3)
+      write (iout,*) "velocity of the center of the mass:"
+      write (iout,*) (vcm(j),j=1,3)
       do j=1,3
         d_t(j,0)=d_t(j,0)-vcm(j)
       enddo
@@ -3028,6 +3073,10 @@
             if (me.eq.king.or..not.out1file) write(iout,*)&
              'Minimizing initial PDB structure: Calling MINIM_DC'
             call minim_dc(etot,iretcode,nfun)
+#ifdef LBFGS
+            if (me.eq.king.or..not.out1file)&
+            write(iout,*) 'LBFGS return code is ',statusbf,' eval ',nfun
+#endif
           else
             call geom_to_var(nvar,varia)
             if(me.eq.king.or..not.out1file) write (iout,*)&
@@ -3036,9 +3085,9 @@
             call var_to_geom(nvar,varia)
 #ifdef LBFGS
             if (me.eq.king.or..not.out1file)&
-            write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
+            write(iout,*) 'LBFGS return code is ',statusbf,' eval ',nfun
             if(me.eq.king.or..not.out1file)&
-            write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
+            write(iout,*) 'LBFGS return code is ',statusbf,' eval ',nfun
 #else
             if (me.eq.king.or..not.out1file)&
             write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
@@ -3240,7 +3289,7 @@
 !      include 'COMMON.TIME1'
       real(kind=8) :: xv,sigv,lowb,highb  ,Ek1
 #ifdef FIVEDIAG
-      integer ichain,n,innt,inct,ibeg,ierr
+      integer ichain,n,innt,inct,ibeg,ierr,innt_org
       real(kind=8) ,allocatable, dimension(:)::  work
       integer,allocatable,dimension(:) :: iwork
 !      double precision Ghalf(mmaxres2_chain),Geigen(maxres2_chain),&
@@ -3253,7 +3302,7 @@
 !#define DEBUG
 #ifdef FIVEDIAG
        real(kind=8) ,allocatable, dimension(:)  :: xsolv,DML,rs
-       real(kind=8) :: sumx,Ek2,Ek3,aux
+       real(kind=8) :: sumx,Ek2,Ek3,aux,masinv
 #ifdef DEBUG
        real(kind=8) ,allocatable, dimension(:)  :: rsold
        real (kind=8),allocatable,dimension(:,:) :: matold,inertia
@@ -3286,15 +3335,21 @@
       EK=0.0d0
       Ek3=0.0d0
 #ifdef DEBUG
-      write(iout,*), nchain
+      write(iout,*), "nchain",nchain
 #endif
       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
+!        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)
+!         innt_org=
+        innt_org=chain_border(1,ichain)
+        if ((molnum(innt_org).eq.5).or.(molnum(innt_org).eq.4)) go to 137
+        if(.not.allocated(ghalf)) print *,"COCO"
+        if(.not.allocated(Ghalf)) allocate(Ghalf(1300*(1300+1)/2))
+        ghalf=0.0d0
         inct=innt+n-1
 #ifdef DEBUG
         write (iout,*) "Chain",ichain," n",n," start",innt
@@ -3364,20 +3419,39 @@
           enddo
         enddo
 #endif
+137     continue
+        write(iout,*) "HERE,",n,innt
+        innt_org=chain_border(1,ichain)
         xv=0.0d0
         ii=0
         do i=1,n
           do k=1,3
             ii=ii+1
+             mnum=molnum(innt_org)
+            if (molnum(innt_org).ge.4) geigen(i)=3.0/msc(itype(innt_org+i-1,mnum),mnum)
+!            if (molnum(innt).eq.5) write(iout,*) "typ",i,innt-1+i,itype(innt+i-1,5)
             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)
+            write (iout,*) "i",i," ii",ii," geigen",geigen(i), &
+           " d_t_work_new",d_t_work_new(ii),innt_org+i-1
           enddo
         enddo
+        if (molnum(innt_org).ge.4) then
+        mnum=molnum(innt_org)
+        do k=1,3
+          do i=1,n
+            ind=(i-1)*3+k
+            d_t_work(ind)=0.0d0
+            masinv=1.0d0/msc(itype(innt_org+i-1,mnum),mnum)
+            d_t_work(ind)=d_t_work(ind)&
+            +masinv*d_t_work_new((i-1)*3+k)
+          enddo
+        enddo
+
+        else
         do k=1,3
           do i=1,n
             ind=(i-1)*3+k
@@ -3386,10 +3460,9 @@
               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
+        endif
 #ifdef DEBUG
         aux=0.0d0
         do k=1,3
@@ -3458,10 +3531,23 @@
         d_t(:,0)=d_t(:,1)
         d_t(:,1)=0.0d0
       endif
+      if (large) then
+        write (iout,*)
+        write (iout,*) "Random vel after 1st transf the Calpha,SC space"
+        write (iout,'(3hORG,1h(,i5,1h),3f10.5)') 0,(d_t(j,0),j=1,3)
+        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
+
 !c      d_a(:,0)=d_a(:,1)
 !c      d_a(:,1)=0.0d0
 !c      write (iout,*) "Shifting accelerations"
       do ichain=2,nchain
+        write(iout,*) "nchain",ichain,chain_border1(1,ichain),molnum(chain_border1(1,ichain))
+        if (molnum(chain_border1(1,ichain)+1).eq.5) cycle
 !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))
@@ -3469,6 +3555,7 @@
       enddo
 !c      write (iout,*) "Adding accelerations"
       do ichain=2,nchain
+        if (molnum(chain_border1(1,ichain)+1).eq.5) cycle
 !c        write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
 !c     &   chain_border(2,ichain-1)
         d_t(:,chain_border1(1,ichain)-1)=&
@@ -3478,10 +3565,23 @@
       do ichain=2,nchain
         write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,&
         chain_border(2,ichain-1)
+        if (molnum(chain_border1(1,ichain)+1).eq.5) cycle
+
         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
+      if (large) then
+        write (iout,*)
+        write (iout,*) "Random vel after 2nd transf the Calpha,SC space"
+        write (iout,'(3hORG,1h(,i5,1h),3f10.5)') 0,(d_t(j,0),j=1,3)
+        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
+
 #else
       ibeg=0
 !c      do j=1,3
@@ -4565,7 +4665,7 @@
        summas=0.0d0
        do i=nnt,nct
          mnum=molnum(i)
-         if ((mnum.ge.5).or.(mnum.eq.3))&
+         if ((mnum.gt.5).or.(mnum.eq.3))&
          mp(mnum)=msc(itype(i,mnum),mnum)
          if (i.lt.nct) then
            summas=summas+mp(mnum)
@@ -4573,11 +4673,11 @@
              vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i))
            enddo
          endif
-         if (mnum.ne.4) then
+!         if (mnum.ne.4) then
          amas=msc(iabs(itype(i,mnum)),mnum)
-         else
-         amas=0.0d0
-         endif
+!         else
+!         amas=0.0d0
+!         endif
 !         amas=msc(iabs(itype(i)))
          summas=summas+amas             
 !         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
@@ -4595,7 +4695,7 @@
            vv(j)=vv(j)+d_t(j,i)
          enddo
        enddo
-!c       write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
+       write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
        do j=1,3
          vcm(j)=vcm(j)/summas
        enddo
@@ -4645,7 +4745,7 @@
         M_PEP=0.0d0
         do i=nnt,nct-1
           mnum=molnum(i)
-          if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
+!          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)
@@ -4660,10 +4760,11 @@
         M_SC=0.0d0		   	 	
         do i=nnt,nct
            mnum=molnum(i)
-           if (mnum.ge.5) cycle
+!           if (mnum.ge.5) cycle
            iti=iabs(itype(i,mnum))	   	 
       	   M_SC=M_SC+msc(iabs(iti),mnum)
            inres=i+nres
+           if (mnum.ge.4) inres=i
            do j=1,3
             cm(j)=cm(j)+msc(iabs(iti),mnum)*c(j,inres)	    
            enddo
@@ -4674,7 +4775,7 @@
 !        write(iout,*) "Center of mass:",cm
         do i=nnt,nct-1
            mnum=molnum(i)
-          if (mnum.ge.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
@@ -4694,8 +4795,9 @@
       	do i=nnt,nct	
            mnum=molnum(i)
            iti=iabs(itype(i,mnum))
-           if (mnum.ge.5) cycle
+!          if (mnum.ge.5) cycle
            inres=i+nres
+           if (mnum.ge.4) inres=i
            do j=1,3
              pr(j)=c(j,inres)-cm(j)	    
            enddo
@@ -4976,7 +5078,7 @@
        summas=0.0d0
        do i=nnt,nct
          mnum=molnum(i)
-         if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
+!         if (mnum.ge.4) mp(mnum)=msc(itype(i,mnum),mnum)
          if (i.lt.nct) then
            summas=summas+mp(mnum)
            do j=1,3
@@ -4984,11 +5086,11 @@
 !             print *,i,j,vv(j),d_t(j,i)
            enddo
          endif
-         if (mnum.ne.4) then 
+!         if (mnum.ne.4) then 
          amas=msc(iabs(itype(i,mnum)),mnum)
-         else
-         amas=0.0d0
-         endif
+!         else
+!         amas=0.0d0
+!         endif
          summas=summas+amas    	     	 
          if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
           .and.(mnum.lt.4)) then
@@ -5004,7 +5106,7 @@
            vv(j)=vv(j)+d_t(j,i)
          enddo
        enddo 
-!       write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
+       write (iout,*) "vcm",(vcm(j),j=1,3)," summas",summas
        do j=1,3
          vcm(j)=vcm(j)/summas
        enddo
@@ -5887,6 +5989,7 @@
 !c          innt=chain_border(1,1)
 !c          inct=chain_border(2,1)
           do i=innt,inct
+            mnum=molnum(i)
             vvec(ind+1)=v_work(j,i)
             ind=ind+1
 !            if (iabs(itype(i)).ne.10) then
@@ -6107,7 +6210,7 @@
 !      include 'COMMON.IOUNITS'
       integer :: IERROR
       integer :: i,j,ind,ind1,m,ichain,innt,inct
-      logical :: lprn = .false.
+      logical :: lprn = .true.
       real(kind=8) :: dtdi !el ,gamvec(2*nres)
 !el      real(kind=8),dimension(2*nres,2*nres) :: ginvfric,fcopy
 !      real(kind=8),allocatable,dimension(:,:) :: fcopy
@@ -6216,12 +6319,11 @@
       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
+        mnum=molnum(i)
           if (iabs(itype(i,1)).ne.10.and.mnum.le.2) then
             ind=ind+2
           else
@@ -6232,11 +6334,11 @@
       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
+         mnum=molnum(i)
           if (iabs(itype(i,1)).ne.10.and.mnum.le.2) then
             DU2fric(ind)=gamvec(i-nnt+1)/4
             DU2fric(ind+1)=0.0d0
@@ -7375,7 +7477,7 @@
       allocate(DM(nres2),DU1(nres2),DU2(nres2))
       allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2))
 !#ifdef DEBUG
-      allocate(Gvec(nres2,nres2))
+      allocate(Gvec(1300,1300))!maximum number of protein in one chain
 !#endif
 #else
       write (iout,*) "Before A Allocation",nfgtasks-1