NEWCORR5D working with 15k, work and iwork in random_vel might need testing
authorAdam Sieradzan <adasko@piasek4.chem.univ.gda.pl>
Mon, 3 Apr 2023 12:25:10 +0000 (14:25 +0200)
committerAdam Sieradzan <adasko@piasek4.chem.univ.gda.pl>
Mon, 3 Apr 2023 12:25:10 +0000 (14:25 +0200)
13 files changed:
PARAM/tube_carbonano2.parm
source/unres/CMakeLists.txt
source/unres/MD.F90
source/unres/REMD.F90
source/unres/compare.F90
source/unres/control.F90
source/unres/energy.F90
source/unres/geometry.F90
source/unres/io.F90
source/unres/io_base.F90
source/unres/io_config.F90
source/unres/minim.F90
source/unres/unres.F90

index ea98cb1..c3e20da 100644 (file)
 1.6606077 5.9355397 0 0 0 0 3.0000000
 1.3228511 5.4343948 0 0 0 0 3.0000000
 1.3228511 5.4343948 0 0 0 0 3.0000000
+   0.000   0.000     0.000     0.000     0.000   0.000E+00  0.00 ! Pro
+   0.000   0.000     0.000     0.000     0.000   0.000E+00  0.00 ! Pro
+   0.000   0.000     0.000     0.000     0.000   0.000E+00  0.00 ! Pro
+   0.000   0.000     0.000     0.000     0.000   0.000E+00  0.00 ! Pro
+   0.000   0.000     0.000     0.000     0.000   0.000E+00  0.00 ! Pro
+   0.000   0.000     0.000     0.000     0.000   0.000E+00  0.00 ! Pro
+   0.000   0.000     0.000     0.000     0.000   0.000E+00  0.00 ! Pro
+
index feeb945..3fa4bc8 100644 (file)
@@ -72,7 +72,7 @@ if (Fortran_COMPILER_NAME STREQUAL "ifort")
   set (CMAKE_Fortran_FLAGS_RELEASE " ")
   set (CMAKE_Fortran_FLAGS_DEBUG   "-O0 -g -traceback")
 #  set(FFLAGS0 "-fpp -c -CB -g -ip " ) 
-  set(FFLAGS0 "-CB -g -ip -fpp  -heap-arrays -mcmodel=large" ) 
+  set(FFLAGS0 "-O3 -g -ip -fpp   -mcmodel=large" ) 
 #  set(FFLAGS0 "-O0 -CB -CA -g" )
   set(FFLAGS1 "-fpp -c -O " ) 
   set(FFLAGS2 "-fpp -c -g -CA -CB ")
index 5a6252e..87f7394 100644 (file)
       do j=1,3
         incr(j)=d_t(j,0)
       enddo
-      do i=nnt,nct-1
+      do i=nnt,nct-1 !czy na pewno nct-1??
        mnum=molnum(i)
 !c        write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3
 !c Skip dummy peptide groups
           do j=1,3
             v(j)=incr(j)+0.5d0*d_t(j,i)
           enddo
-          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,*) "Kinetic trp:",i,(v(j),j=1,3)
           vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
           KEt_p=KEt_p+mp(mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
       do i=nnt,nct
        mnum=molnum(i)
         iti=iabs(itype(i,mnum))
+        if (mnum.eq.5) iti=itype(i,mnum)
         if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
            .or.mnum.ge.3) then
           do j=1,3
             v(j)=incr(j)+d_t(j,nres+i)
          enddo
         endif
-        write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
-        write (iout,*) "i",i," msc",msc(iti,1)," v",(v(j),j=1,3)
+!        if (mnum.ne.5) then
+!        write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
+!        write (iout,*) "i",i," msc",msc(iti,mnum)," 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))
         vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
+!        endif
         do j=1,3
           incr(j)=incr(j)+d_t(j,i)
         enddo
        enddo
 !c The total kinetic energy      
   111  continue
-!c       write(iout,*) ' KEt_p',KEt_p,' KEt_sc',KEt_sc,' KEr_p',KEr_p,
-!c     &  ' KEr_sc', KEr_sc
+!       write(iout,*) ' KEt_p',KEt_p,' KEt_sc',KEt_sc,' KEr_p',KEr_p, &
+!       ' KEr_sc', KEr_sc
        KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc)
 !c       write (iout,*) "KE_total",KE_total
        return
       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)
 !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))
         endif
         if (mod(itime,ntwx).eq.0) then
           call returnbox
+          call enerprint(potEcomp)
+
           write (tytul,'("time",f8.2)') totT
           if(mdpdb) then
              call hairpin(.true.,nharp,iharp)
              call secondary2(.true.)
              call pdbout(potE,tytul,ipdb)
-             call enerprint(potEcomp)
+!             call enerprint(potEcomp)
           else 
              call cartout(totT)
           endif
       if(me.eq.king.or..not.out1file)then
        write (iout,*) "Initial velocities"
        do i=0,nres
-         write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
+         write (iout,'(i6,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
          (d_t(j,i+nres),j=1,3)
        enddo
 !  Zeroing the total angular momentum of the system
         write (iout,*) "Initial velocities"
         write (iout,"(13x,' backbone ',23x,' side chain')")
         do i=0,nres
-          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
+          write (iout,'(i6,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
           (d_t(j,i+nres),j=1,3)
         enddo
         write (iout,*) "Initial accelerations"
       real(kind=8) :: xv,sigv,lowb,highb  ,Ek1
 #ifdef FIVEDIAG
       integer ichain,n,innt,inct,ibeg,ierr
-      double precision work(48*nres)
-      integer iwork(6*nres)
+      real(kind=8) ,allocatable, dimension(:)::  work
+      integer,allocatable,dimension(:) :: iwork
 !      double precision Ghalf(mmaxres2_chain),Geigen(maxres2_chain),&
 !      Gvec(maxres2_chain,maxres2_chain)
 !      common /przechowalnia/Ghalf,Geigen,Gvec
-#ifdef DEBUG
-      double precision inertia(maxres2_chain,maxres2_chain)
-#endif
+!#ifdef DEBUG
+!      double precision inertia(maxres2_chain,maxres2_chain)
+!#endif
 #endif
-#define DEBUG
+!#define DEBUG
 #ifdef FIVEDIAG
-       real(kind=8) ,allocatable, dimension(:)  :: DDU1,DDU2,DL2,DL1,xsolv,DML,rs
+       real(kind=8) ,allocatable, dimension(:)  :: xsolv,DML,rs
        real(kind=8) :: sumx,Ek2,Ek3,aux
 #ifdef DEBUG
        real(kind=8) ,allocatable, dimension(:)  :: rsold
        real (kind=8),allocatable,dimension(:,:) :: matold,inertia
        integer :: iti
 #endif
-      integer :: i,j,ii,k,ind,mark,imark,mnum
+#endif
+      integer :: i,j,ii,k,mark,imark,mnum,nres2
+      integer(kind=8) :: ind
 ! Generate random velocities from Gaussian distribution of mean 0 and std of KT/m 
+!#undef DEBUG
 ! First generate velocities in the eigenspace of the G matrix
 !      write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
 !      call flush(iout)
-!#define DEBUG
+#ifdef FIVEDIAG
+       if(.not.allocated(work)) then
+       allocate(work(48*nres))
+       allocate(iwork(6*nres))
+       endif
+       print *,"IN RANDOM VEL"
+       nres2=2*nres
+!       print *,size(ghalf)
+#undef DEBUG
 #ifdef DEBUG
       write (iout,*) "Random_vel, fivediag"
+      flush(iout)
       allocate(inertia(2*nres,2*nres))
 #endif
       d_t=0.0d0
 #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
         n=dimen_chain(ichain)
         innt=iposd_chain(ichain)
         mnum1=molnum(i+1)
         if (itype(i,mnum).eq.ntyp1_molec(mnum)&
          .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
-          if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
+!          if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
+          if (mnum.ge.5) mp(mnum)=0.0d0
           M_PEP=M_PEP+mp(mnum)
 
         do j=1,3
       do i=nnt,nct-1
         mnum=molnum(i)
         mnum1=molnum(i+1)
-        if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
+!        if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
+          if (mnum.ge.5) mp(mnum)=0.0d0
         if (itype(i,mnum).eq.ntyp1_molec(mnum)&
          .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
         do j=1,3
       do i=nnt,nct-1
         mnum=molnum(i)
         mnum1=molnum(i+1)
-        if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
+!        if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
+          if (mnum.ge.5) mp(mnum)=0.0d0
         if (itype(i,mnum).eq.ntyp1_molec(mnum)&
         .or. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) cycle
         do j=1,3
         call vecpr(pr(1),v(1),vp)
 !c          write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),
 !c      &     " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
-        if (mnum.gt.4) then
-         mscab=0.0d0
-        else
+!        if (mnum.gt.4) then
+!         mscab=0.0d0
+!        else
          mscab=msc(iti,mnum)
-        endif
+!        endif
         do j=1,3
           L(j)=L(j)+mscab*vp(j)
         enddo
 ! commom.langevin.lang0
 !      common /langforc/
       allocate(friction(3,0:nres2),stochforc(3,0:nres2)) !(3,0:MAXRES2)
+#ifndef FIVEDIAG
       if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2))
       allocate(fricvec(nres2,nres2)) !(MAXRES2,MAXRES2)
+#endif
       allocate(fric_work(nres6),stoch_work(nres6),fricgam(nres6)) !(MAXRES6)
       allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch)
 #endif
-
+#ifndef FIVEDIAG
       if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2))
+#endif
 !----------------------
 ! commom.hairpin in CSA module
 !----------------------
 #ifdef FIVEDIAG
       allocate(DM(nres2),DU1(nres2),DU2(nres2))
       allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2))
+!#ifdef DEBUG
       allocate(Gvec(nres2,nres2))
+!#endif
 #else
       write (iout,*) "Before A Allocation",nfgtasks-1
       call flush(iout)
index e9df38f..8649b12 100644 (file)
@@ -24,7 +24,9 @@
 !       include 'DIMENSIONS'
       use comm_cipiszcze
       use energy_data
+#ifdef FIVEDIAG
       use energy, only: grad_transform
+#endif
       use geometry_data, only: nres
       use control_data    !el, only: mucadyn,lmuca
 #ifdef MPI
@@ -51,7 +53,7 @@
        real(kind=8) :: rscheck(2*nres),rsold(2*nres)
 #endif
 #endif
-       logical :: lprn = .true.
+       logical :: lprn = .false.
 !el       common /cipiszcze/ itime
        itime = itt_comm
 
         enddo
         inct_prev=inct+1
         do i=innt,inct
-          if (itype(i).ne.10) then
+          mnum=molnum(i)
+          if ((itype(i).ne.10).and.(mnum.le.3)) then
             do j=1,3
               d_a(j,i+nres)=d_a(j,i+nres)-d_a(j,i)
             enddo
       include 'mpif.h'
       integer :: ierror
       real(kind=8) :: time00
-      real(kind=8) ,allocatable, dimension(:)  :: DDM,DDU1,DDU2
+!      real(kind=8) ,allocatable, dimension(:)  :: DDM,DDU1,DDU2
 #endif
 !      include 'COMMON.SETUP'
 !      include 'COMMON.VAR'
       real(kind=8),allocatable,dimension(:,:) :: Bmat,matmult
 #endif
       real(kind=8) :: dtdi
-      real(kind=8),dimension(2*nres) :: massvec,sqreig !(maxres2) maxres2=2*maxres
+      real(kind=8),dimension(:),allocatable :: massvec,sqreig  !(maxres2) maxres2=2*maxres
       real(kind=8) :: relfeh,eps1,eps2
 !el      real(kind=8),dimension(:),allocatable :: Ghalf
 !el      real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf  !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2))
 !el      real(kind=8),dimension(2*nres,2*nres) :: Gcopy        !(maxres2,maxres2)
 !el      real(kind=8),dimension(:,:),allocatable :: Gcopy
-      real(kind=8),dimension(8*6*nres) :: work !(8*maxres6)
-      integer,dimension(6*nres) :: iwork       !(maxres6) maxres6=6*maxres
+      real(kind=8),dimension(:),allocatable :: work    !(8*maxres6)
+      integer,dimension(:),allocatable :: iwork        !(maxres6) maxres6=6*maxres
+!      common /jakistam/ iwork,work,massvec,sqreig
 !el      common /przechowalnia/ Gcopy,Ghalf
       real(kind=8) :: coeff,mscab
       integer :: i,j,ind,ind1,k,l,ii,jj,m,m1,ii1,iti,nres2,ierr,nind,mark
       real(kind=8) :: ip4
       integer :: iz,mnum,ichain,n,dimenp,innt,inct
+      if(.not.allocated(massvec)) then
+      allocate(massvec(2*nres),sqreig(2*nres))
+      allocate(work(8*6*nres))
+      allocate(iwork(6*nres))
+      endif
       print *,"just entered"
       relfeh=1.0d-14
       nres2=2*nres
       print *,"FIVEDIAG"
       write (iout,*) "before FIVEDIAG"
+#ifdef FIVEDIAG
       if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !mmaxres2=(maxres2*(maxres+1)/2)
+#endif
 #ifndef FIVEDIAG
+      if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !mmaxres2=(maxres2*(maxres+1)/2)
       write (iout,*) "ALLOCATE"
       print *,"ALLOCATE"
       if(.not.allocated(Gcopy)) allocate(Gcopy(nres2,nres2)) !(maxres2,maxres2)
         innt=chain_border(1,ichain)
         mnum=molnum(innt)
         inct=chain_border(2,ichain)
+        if (mnum.eq.5) mp(mnum)=0.0
+!        if (mnum.eq.5) mp(mnum)=msc(itype(innt,mnum),mnum)
         DM(ind)=mp(mnum)/4+ip(mnum)/4
         if (iabs(itype(innt,1)).eq.10.or.molnum(innt).gt.2) then
-          DM(ind)=DM(ind)+msc(itype(innt,mnum),mnum)
+           DM(ind)=DM(ind)+msc(itype(innt,mnum),mnum)
           ind=ind+1
           nind=1
         else
         do i=innt+1,inct-1
          mnum=molnum(i)
 !        if (iabs(itype(i)).eq.ntyp1) cycle
-          if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
-          if (mnum.eq.5) ip(mnum)=isc(itype(i,mnum),mnum)
+!          if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
+!          if (mnum.eq.5) ip(mnum)=isc(itype(i,mnum),mnum)
+        if (mnum.eq.5) mp(mnum)=0.0
+!        if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
           DM(ind)=2*ip(mnum)/4+mp(mnum)/2
           if (iabs(itype(i,1)).eq.10.or.molnum(i).gt.2) then
-            if (iabs(itype(i,1)).eq.10.or.molnum(i).gt.2)&
-             DM(ind)=DM(ind)+msc(itype(i,molnum(i)),mnum)
+!            if (iabs(itype(i,1)).eq.10.or.molnum(i).gt.2)&
+!             DM(ind)=DM(ind)+msc(itype(i,molnum(i)),mnum)
+             DM(ind)=DM(ind)+msc(itype(i,mnum),mnum)
             ind=ind+1
             nind=nind+1
           else
         if (inct.gt.innt) then
 !         DM(ind)=ip4+mp(molnum(inct))/4
           mnum=molnum(inct)
+        if (mnum.eq.5) mp(mnum)=0.0
+!        if (mnum.eq.5) mp(mnum)=msc(itype(inct,molnum(inct)),molnum(inct))
+
           DM(ind)=mp(mnum)/4+ip(mnum)/4
           if (iabs(itype(inct,mnum)).eq.10.or.molnum(inct).gt.2) then
-            DM(ind)=DM(ind)+msc(itype(inct,molnum(inct)),molnum(inct))
+              DM(ind)=DM(ind)+msc(itype(inct,molnum(inct)),molnum(inct))
             ind=ind+1
             nind=nind+1
           else
             DU1(ind+1)=0.0d0
             ind=ind+2
           else
-            if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
-            if (mnum.eq.5) ip(mnum)=isc(itype(i,mnum),mnum)
+!           if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
+!            if (mnum.eq.5) ip(mnum)=isc(itype(i,mnum),mnum)
+             if (mnum.eq.5) mp(mnum)=0.0
             DU1(ind)=mp(mnum)/4-ip(mnum)/4
             ind=ind+1
           endif
 !       if (iabs(itype(i)).eq.ntyp1) cycle
 !c          write (iout,*) "i",i," itype",itype(i),ntyp1
           if (iabs(itype(i,1)).ne.10 .and. iabs(itype(i,mnum)).ne.ntyp1_molec(mnum).and.mnum.le.2) then
-            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.0
             DU2(ind)=mp(mnum)/4-ip(mnum)/4
             DU2(ind+1)=0.0d0
             ind=ind+2
index 6ad8870..2f47cf6 100644 (file)
       if(.not.allocated(bfrag)) allocate(bfrag(4,nres/3)) !(4,maxres/3)
       if(.not.allocated(hfrag)) allocate(hfrag(2,nres/3)) !(2,maxres/3)
 !      COMMON /WAGI/
+#ifndef FIVEDIAG
       if(.not.allocated(w)) allocate(w(maxres22),d0(maxres22)) !(maxres22)
 !      COMMON /POCHODNE/
 !el      allocate(DRDG(maxres22,maxres22)) !(MAXRES22,MAXRES)
       if (.not.allocated(DDD)) allocate(DDD(maxres22)) !(maxres22)
       if (.not.allocated(H)) allocate(H(nres,nres)) !(MAXRES,MAXRES)
+#endif
       if (.not.allocated(XX)) allocate(XX(nres)) !(MAXRES)
 !      COMMON /frozen/
       if (.not.allocated(mask)) allocate(mask(nres)) !(maxres)
index d6293ba..516a8d7 100644 (file)
             my_ele_inds_vdw_nucl,my_ele_inde_vdw_nucl,ind_eleint_vdw_nucl,&
             ind_eleint_vdw_old_nucl,nscp_int_tot_nucl,my_scp_inds_nucl,&
             my_scp_inde_nucl,ind_scpint_nucl,ind_scpint_old_nucl,impishi
-       integer,dimension(nres,nres) :: remmat
+       integer(kind=1),dimension(:,:),allocatable :: remmat
 !      integer,dimension(5) :: nct_molec,nnt_molec
 !el      allocate(itask_cont_from(0:nfgtasks-1)) !(0:max_fg_procs-1)
 !el      allocate(itask_cont_to(0:nfgtasks-1)) !(0:max_fg_procs-1)
 !... Determine the numbers of start and end SC-SC interaction
 !... to deal with by current processor.
 !write (iout,*) '******INIT_INT_TABLE nres=',nres,' nnt=',nnt,' nct=',nct
+      print *,"in spliting contacts"
       do i=0,nfgtasks-1
         itask_cont_from(i)=fg_rank
         itask_cont_to(i)=fg_rank
       lprint=energy_dec
       itmp=0
       do i=1,5
+       print *,i,nres_molec(i)
        if (nres_molec(i).eq.0) cycle
       itmp=itmp+nres_molec(i)
       if (itype(itmp,i).eq.ntyp1_molec(i)) then
       nnt_molec(i)=itmp+1
       endif
       enddo
+!      if (.not.allocated(nres_molec)) print *,"WHATS WRONG"
       print *,"nres_molec",nres_molec(:)
       print *,"nnt_molec",nnt_molec(:)
       print *,"nct_molec",nct_molec(:)
        ilipang_end=itmp+nres_molec(4)-1
 !      create LJ LIST MAXIMUM
 !      Eliminate branching from list
+       if(.not.allocated(remmat))&
+        allocate(remmat(itmp+1:nres_molec(4)+itmp,itmp+1:nres_molec(4)+itmp))
           remmat=0
        do i=1+itmp,nres_molec(4)-1+itmp
         if (itype(i,4).eq.12) ibra=i
 !        write (iout,*) "Loop backward iatel_e-1",iatel_e-1,
 !     &     " iatel_e",iatel_e
 !        call flush(iout)
+#ifndef NEWCORR
         nat_sent=0
         do i=iatel_s,iatel_e
 !          write (iout,*) "i",i," ielstart",ielstart(i),
             iat_sent(nat_sent)=i
           endif
         enddo
+#endif
         if (lprint) then
         write (iout,*)"After longrange ntask_cont_from",ntask_cont_from,&
          " ntask_cont_to",ntask_cont_to
         write (iout,*) "itask_cont_to",&
           (itask_cont_to(i),i=1,ntask_cont_to)
         call flush(iout)
+#ifndef NEWCORR
         write (iout,*) "iint_sent"
         do i=1,nat_sent
           ii=iat_sent(i)
           write (iout,'(20i4)') ii,(j,(iint_sent(k,j,i),k=1,4),&
             j=ielstart(ii),ielend(ii))
         enddo
+#endif
         write (iout,*) "iturn3_sent iturn3_start",iturn3_start,&
           " iturn3_end",iturn3_end
         write (iout,'(20i4)') (i,(iturn3_sent(j,i),j=1,4),&
           itask_cont_from(0),CONT_FROM_GROUP,IERR)
         call MPI_Group_incl(fg_group,ntask_cont_to+1,itask_cont_to(0),&
           CONT_TO_GROUP,IERR)
+#ifndef NEWCORR
         do i=1,nat_sent
           ii=iat_sent(i)
           iaux=4*(ielend(ii)-ielstart(ii)+1)
 !          write (iout,*) "Ranks translated i=",i
 !          call flush(iout)
         enddo
+#endif
         iaux=4*(iturn3_end-iturn3_start+1)
           if (iaux.lt.0) iaux=0
         call MPI_Group_translate_ranks(fg_group,iaux,&
            iturn4_sent(1,iturn4_start),CONT_TO_GROUP,&
            iturn4_sent_local(1,iturn4_start),IERR)
         if (lprint) then
+#ifndef NEWCORR
+
         write (iout,*) "iint_sent_local"
         do i=1,nat_sent
           ii=iat_sent(i)
             j=ielstart(ii),ielend(ii))
           call flush(iout)
         enddo
+#endif
         if (iturn3_end.gt.0) then
         write (iout,*) "iturn3_sent_local iturn3_start",iturn3_start,&
           " iturn3_end",iturn3_end
index 4dc58b1..d81310f 100644 (file)
 !-----------------------------------------------------------------------------
 ! common.sbridge
 !      common /dyn_ssbond/
-      real(kind=8),dimension(:,:),allocatable :: dyn_ssbond_ij !(maxres,maxres)
+      real(kind=8),dimension(:),allocatable :: dyn_ssbond_ij !(maxres,maxres)
 !-----------------------------------------------------------------------------
 ! common.sccor
 ! Parameters of the SCCOR term
 #ifdef MPI      
       real(kind=8) :: weights_(n_ene) !,time_Bcast,time_Bcastw
 ! shielding effect varibles for MPI
-      real(kind=8) ::  fac_shieldbuf(nres), &
-      grad_shield_locbuf1(3*maxcontsshi*nres), &
-      grad_shield_sidebuf1(3*maxcontsshi*nres), &
-      grad_shield_locbuf2(3*maxcontsshi*nres), &
-      grad_shield_sidebuf2(3*maxcontsshi*nres), &
-      grad_shieldbuf1(3*nres), &
-      grad_shieldbuf2(3*nres)
-
-       integer ishield_listbuf(-1:nres), &
-       shield_listbuf(maxcontsshi,-1:nres),k,j,i,iii,impishi,mojint,jjj
+      real(kind=8) ::  fac_shieldbuf(nres_molec(1)), &
+      grad_shield_locbuf1(3*maxcontsshi*nres_molec(1)), &
+      grad_shield_sidebuf1(3*maxcontsshi*nres_molec(1)), &
+      grad_shield_locbuf2(3*maxcontsshi*nres_molec(1)), &
+      grad_shield_sidebuf2(3*maxcontsshi*nres_molec(1)), &
+      grad_shieldbuf1(3*nres_molec(1)), &
+      grad_shieldbuf2(3*nres_molec(1))
+
+       integer ishield_listbuf(-1:nres_molec(1)), &
+       shield_listbuf(maxcontsshi,-1:nres_molec(1)),k,j,i,iii,impishi,mojint,jjj
 !       print *,"I START ENERGY"
        imatupdate=100
 !       if (mod(itime_mat,imatupdate).eq.0) call make_SCSC_inter_list
         time_Bcastw=time_Bcastw+MPI_Wtime()-time00
 !        call chainbuild_cart
       endif
-!       print *,"itime_mat",itime_mat,imatupdate
+       print *,"itime_mat",itime_mat,imatupdate
         if (nfgtasks.gt.1) then 
         call MPI_Bcast(itime_mat,1,MPI_INT,king,FG_COMM,IERROR)
         endif
         Eafmforce=0.0d0
       endif
       endif
-!      print *,"before tubemode",tubemode
+      print *,"before tubemode",tubemode
       if (tubemode.eq.1) then
        call calctube(etube)
       else if (tubemode.eq.2) then
       else
        etube=0.0d0
       endif
-!      print *, "TU JEST PRZED EHPB"
+      print *, "TU JEST PRZED EHPB"
       call edis(ehpb)
 
 !--------------------------------------------------------
 !      include 'COMMON.SBRIDGE'
       logical :: lprn
 !el local variables
-      integer :: iint,itypi,itypi1,itypj,subchap,icont
+      integer :: iint,itypi,itypi1,itypj,subchap,icont,countss
       real(kind=8) :: rrij,xi,yi,zi,sig,rij_shift,fac,e1,e2,sigm,epsi
       real(kind=8) :: evdw,sig0ij
       real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
 !     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
       lprn=.false.
+      countss=0
 !     if (icall.eq.0) lprn=.false.
 !el      ind=0
       dCAVdOM2=0.0d0
 !        do iint=1,nint_gr(i)
 !          do j=istart(i,iint),iend(i,iint)
             IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
-              call dyn_ssbond_ene(i,j,evdwij)
+              countss=countss+1
+              call dyn_ssbond_ene(i,j,evdwij,countss)
               evdw=evdw+evdwij
               if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') &
                               'evdw',i,j,evdwij,' ss'
@@ -15207,7 +15209,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!!
 !      include 'COMMON.CONTROL'
       logical :: lprn
 !el local variables
-      integer :: iint,itypi,itypi1,itypj,subchap
+      integer :: iint,itypi,itypi1,itypj,subchap,countss
       real(kind=8) :: rrij,xi,yi,zi,fac,sigm,epsi,sig0ij,sig
       real(kind=8) :: sss,e1,e2,evdw,rij_shift,sss_grad
       real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
@@ -15218,6 +15220,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!!
 !     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
       lprn=.false.
+      countss=0
 !     if (icall.eq.0) lprn=.false.
 !el      ind=0
       do i=iatsc_s,iatsc_e
@@ -15244,7 +15247,8 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!!
         do iint=1,nint_gr(i)
           do j=istart(i,iint),iend(i,iint)
             IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
-              call dyn_ssbond_ene(i,j,evdwij)
+              countss=countss+1
+              call dyn_ssbond_ene(i,j,evdwij,countss)
               evdw=evdw+evdwij
               if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') &
                               'evdw',i,j,evdwij,' ss'
@@ -19571,14 +19575,14 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!!
 !EL      external ran_number
 
 !     Local variables
-      integer :: i,j,k,l,lmax,p,pmax
+      integer :: i,j,k,l,lmax,p,pmax,countss
       real(kind=8) :: rmin,rmax
       real(kind=8) :: eij
 
       real(kind=8) :: d
       real(kind=8) :: wi,rij,tj,pj
 !      return
-
+      countss=1
       i=5
       j=14
 
@@ -19625,14 +19629,14 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!!
             dc_norm(k,nres+j)=dc(k,nres+j)/d
          enddo
 
-         call dyn_ssbond_ene(i,j,eij)
+         call dyn_ssbond_ene(i,j,eij,countss)
       enddo
       enddo
       call exit(1)
       return
       end subroutine check_energies
 !-----------------------------------------------------------------------------
-      subroutine dyn_ssbond_ene(resi,resj,eij)
+      subroutine dyn_ssbond_ene(resi,resj,eij,countss)
 !      implicit none
 !      Includes
       use calc_data
@@ -19665,7 +19669,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!!
 
 !     Local variables
       logical :: havebond
-      integer itypi,itypj
+      integer itypi,itypj,countss
       real(kind=8) :: rrij,ssd,deltat1,deltat2,deltat12,cosphi
       real(kind=8) :: sig0ij,ljd,sig,fac,e1,e2
       real(kind=8),dimension(3) :: dcosom1,dcosom2
@@ -19963,9 +19967,9 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!!
 !        endif
 !#endif
 !#endif
-      dyn_ssbond_ij(i,j)=eij
-      else if (.not.havebond .and. dyn_ssbond_ij(i,j).lt.1.0d300) then
-      dyn_ssbond_ij(i,j)=1.0d300
+      dyn_ssbond_ij(countss)=eij
+      else if (.not.havebond .and. dyn_ssbond_ij(countss).lt.1.0d300) then
+      dyn_ssbond_ij(countss)=1.0d300
 !#ifndef CLUST
 !#ifndef WHAM
 !        write(iout,'(a15,f12.2,f8.1,2i5)')
@@ -20246,10 +20250,10 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!!
 !      include 'COMMON.MD'
 !     Local variables
       real(kind=8) :: emin
-      integer :: i,j,imin,ierr
+      integer :: i,j,imin,ierr,k
       integer :: diff,allnss,newnss
       integer,dimension(maxdim) :: allflag,allihpb,alljhpb,& !(maxdim)(maxdim=(maxres-1)*(maxres-2)/2)
-            newihpb,newjhpb
+            newihpb,newjhpb,aliass
       logical :: found
       integer,dimension(0:nfgtasks) :: i_newnss
       integer,dimension(0:nfgtasks) :: displ
@@ -20257,14 +20261,19 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!!
       integer :: g_newnss
 
       allnss=0
+      k=0
       do i=1,nres-1
       do j=i+1,nres
-        if (dyn_ssbond_ij(i,j).lt.1.0d300) then
+        if ((itype(i,1).eq.1).and.(itype(j,1).eq.1)) then
+        k=k+1
+        if (dyn_ssbond_ij(k).lt.1.0d300) then
           allnss=allnss+1
           allflag(allnss)=0
           allihpb(allnss)=i
           alljhpb(allnss)=j
-        endif
+          aliass(allnss)=k
+       endif
+       endif
       enddo
       enddo
 
@@ -20273,8 +20282,8 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!!
  1    emin=1.0d300
       do i=1,allnss
       if (allflag(i).eq.0 .and. &
-           dyn_ssbond_ij(allihpb(i),alljhpb(i)).lt.emin) then
-        emin=dyn_ssbond_ij(allihpb(i),alljhpb(i))
+           dyn_ssbond_ij(aliass(allnss)).lt.emin) then
+        emin=dyn_ssbond_ij(aliass(allnss))
         imin=i
       endif
       enddo
@@ -20976,7 +20985,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!!
 !C         fac=fac+faccav
 !C 667     continue
        endif
-        if (energy_dec) write(iout,*),i,rdiff,enetube(i),enecavtube(i)
+        if (energy_dec) write(iout,*),"ETUBE_PEP",i,rdiff,enetube(i),enecavtube(i)
       do j=1,3
       gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac/2.0d0
       gg_tube(j,i)=gg_tube(j,i)+vectube(j)*fac/2.0d0
@@ -21052,7 +21061,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!!
         gg_tube_SC(j,i)=gg_tube_SC(j,i)+vectube(j)*fac
         gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac
        enddo
-        if (energy_dec) write(iout,*),i,rdiff,enetube(i+nres),enecavtube(i+nres)
+        if (energy_dec) write(iout,*),"ETUBE",i,rdiff,enetube(i+nres),enecavtube(i+nres)
       enddo
 
       
@@ -21109,7 +21118,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!!
        do j=1,3
         gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac
        enddo
-        if (energy_dec) write(iout,*),i,rdiff,enetube(i+nres)
+        if (energy_dec) write(iout,*) "ETUBLIP",i,rdiff,enetube(i+nres)
       enddo           
 
 
@@ -21547,6 +21556,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!!
 !      common /contacts1/
       allocate(num_cont(0:nres+4))
 !(maxres)
+#ifndef NEWCORR
       allocate(jcont(maxconts,nres))
 !(maxconts,maxres)
       allocate(facont(maxconts,nres))
@@ -21571,9 +21581,10 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!!
       allocate(ees0plist(maxconts,nres))
       
 !(maxconts,maxres)
-      allocate(num_cont_hb(nres))
 !(maxres)
       allocate(jcont_hb(maxconts,nres))
+#endif
+      allocate(num_cont_hb(nres))
 !(maxconts,maxres)
 !      common /rotat/
       allocate(Ug(2,2,nres))
@@ -21648,9 +21659,9 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!!
       allocate(sintab2(nres))
 !(maxres)
 !      common /dipmat/ 
-      allocate(a_chuj(2,2,maxconts,nres))
+!      allocate(a_chuj(2,2,maxconts,nres))
 !(2,2,maxconts,maxres)(maxconts=maxres/4)
-      allocate(a_chuj_der(2,2,3,5,maxconts,nres))
+!      allocate(a_chuj_der(2,2,3,5,maxconts,nres))
 !(2,2,3,5,maxconts,maxres)(maxconts=maxres/4)
 !      common /contdistrib/
       allocate(ncont_sent(nres))
@@ -21658,8 +21669,12 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!!
 
       allocate(iat_sent(nres))
 !(maxres)
+#ifndef NEWCORR
+      print *,"before iint_sent allocate"
       allocate(iint_sent(4,nres,nres))
       allocate(iint_sent_local(4,nres,nres))
+      print *,"after iint_sent allocate"
+#endif
 !(4,maxres,maxres)
       allocate(iturn3_sent(4,0:nres+4))
       allocate(iturn4_sent(4,0:nres+4))
@@ -21675,8 +21690,15 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!!
 !----------------------
 ! commom.deriv;
 !      common /derivat/ 
+#ifdef NEWCORR
+      print *,"before dcdv allocate"
+      allocate(dcdv(6,nres+2))
+      allocate(dxdv(6,nres+2))
+#else
+      print *,"before dcdv allocate"
       allocate(dcdv(6,maxdim))
       allocate(dxdv(6,maxdim))
+#endif
 !(6,maxdim)
       allocate(dxds(6,nres))
 !(6,maxres)
@@ -21871,11 +21893,11 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!!
 !el      integer,dimension(:),allocatable :: ihpb,jhpb,ibecarb !(maxdim) !el ibecarb !!! nie używane
 !      common /dyn_ssbond/
 ! and side-chain vectors in theta or phi.
-      allocate(dyn_ssbond_ij(0:nres+4,0:nres+4))
+      allocate(dyn_ssbond_ij(10000))
 !(maxres,maxres)
 !      do i=1,nres
 !        do j=i+1,nres
-      dyn_ssbond_ij(:,:)=1.0d300
+      dyn_ssbond_ij(:)=1.0d300
 !        enddo
 !      enddo
 
@@ -21934,6 +21956,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!!
       allocate(uygrad(3,3,2,nres))
       allocate(uzgrad(3,3,2,nres))
 !(3,3,2,maxres)
+      print *,"before all 300"
 ! allocateion of lists JPRDLA
       allocate(newcontlistppi(300*nres))
       allocate(newcontlistscpi(350*nres))
@@ -26751,7 +26774,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!!
 !      include 'COMMON.SBRIDGE'
       logical :: lprn
 !el local variables
-      integer :: iint,itypi1,subchap,isel
+      integer :: iint,itypi1,subchap,isel,countss
       real(kind=8) :: rrij,xi,yi,zi,sig,rij_shift,e1,e2,sigm,epsi
       real(kind=8) :: evdw,aa,bb
       real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
@@ -26776,6 +26799,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!!
        evdw=0.0d0
        eps_out=80.0d0
        sss_ele_cut=1.0d0
+       countss=0
 !       print *,"EVDW KURW",evdw,nres
       do i=iatsc_s,iatsc_e
 !        print *,"I am in EVDW",i
@@ -26804,7 +26828,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!!
         do j=istart(i,iint),iend(i,iint)
 !             print *,"JA PIER",i,j,iint,istart(i,iint),iend(i,iint)
           IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
-            call dyn_ssbond_ene(i,j,evdwij)
+            call dyn_ssbond_ene(i,j,evdwij,countss)
             evdw=evdw+evdwij
             if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') &
                         'evdw',i,j,evdwij,' ss'
index 60c888f..866f5e2 100644 (file)
         write (iout,1212) restyp(itype(i,1),1),i,dist(i-1,i),&
         alfai,be,dist(nres+i,i),rad2deg*alpha(nres+i,i,nres2+2),be1
       enddo   
- 1212 format (a3,'(',i3,')',2(f10.5,2f10.2))
+ 1212 format (a3,'(',i6,')',2(f10.5,2f10.2))
 
       endif
 
        rad2deg*alph(i),rad2deg*omeg(i)
       enddo
       endif
- 1212 format (a3,'(',i3,')',2(f15.10,2f10.2))
+ 1212 format (a3,'(',i6,')',2(f15.10,2f10.2))
 #ifdef TIMING
       time_intfcart=time_intfcart+MPI_Wtime()-time01
 #endif
                 '     centroid coordinates'/ &
                 '       ', 6X,'X',11X,'Y',11X,'Z',&
                                 10X,'X',11X,'Y',11X,'Z')
-  110 format (a,'(',i3,')',6f12.5)
+  110 format (a,'(',i6,')',6f12.5)
       return
       end subroutine cartprint
 !-----------------------------------------------------------------------------
 !      common /rotmat/
       allocate(t(3,3,nres),r(3,3,nres))
       allocate(prod(3,3,nres),rt(3,3,nres)) !(3,3,maxres)
+      print *,"before permutations",maxperm,maxsym
 !      common /refstruct/
       if(.not.allocated(cref)) allocate(cref(3,nres2+2,maxperm)) !(3,maxres2+2,maxperm)
+      print *,"cref"
 !elwrite(iout,*) "jestem w alloc geo 2"
 !      allocate(crefjlee(3,nres2+2)) !(3,maxres2+2)
         if (.not.allocated(crefjlee)) allocate (crefjlee(3,nres2+2))
       if(.not.allocated(chain_rep)) allocate(chain_rep(3,nres2+2,maxsym)) !(3,maxres2+2,maxsym)
+      print *,"chain_rep"
       if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
 !      common /from_zscore/ in module.compare
 !----------------------
index 8ec3d19..4750e57 100644 (file)
 
 !        print *,'Finished reading pdb data'
         if(me.eq.king.or..not.out1file) &
-         write (iout,'(a,i3,a,i3)')'nsup=',nsup,&
+         write (iout,'(a,i7,a,i7)')'nsup=',nsup,&
          ' nstart_sup=',nstart_sup !,"ergwergewrgae"
 !el        if(.not.allocated(itype_pdb)) 
         allocate(itype_pdb(nres))
 !       print *,"tu dochodze"
 ! znamy nres oraz nss można zaalokowac potrzebne tablice
       call alloc_geo_arrays
+      print *,"after  alloc_geo_arrays"
       call alloc_ener_arrays
+      print *,"after  alloc_ener_arrays"
 !--------------------------------
 ! 8/13/98 Set limits to generating the dihedral angles
       do i=1,nres
       if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1
       if (pdbref) then
         if(me.eq.king.or..not.out1file) &
-         write (iout,'(a,i3)') 'nsup=',nsup
+         write (iout,'(a,i7)') 'nsup=',nsup
         nstart_seq=nnt
         if (nsup.le.(nct-nnt+1)) then
           do i=0,nct-nnt+1-nsup
index 21ecd0e..1d7c037 100644 (file)
@@ -5,7 +5,7 @@
       implicit none
 !-----------------------------------------------------------------------------
 ! Max. number of AA residues
-      integer,parameter :: maxres=6500!1200
+      integer,parameter :: maxres=101000!1200
 ! Appr. max. number of interaction sites
       integer,parameter :: maxres2=2*maxres
 !      parameter (maxres6=6*maxres)
       do i=1,nres
        iti=itype(i,1)
 !         print *,vbld(i),"vbld(i)",i
-        write (iout,'(a3,i4,6f10.3)') restyp(iti,1),i,vbld(i),&
+        write (iout,'(a3,i6,6f10.3)') restyp(iti,1),i,vbld(i),&
            rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),rad2deg*alph(i),&
            rad2deg*omeg(i)
       enddo
index 4d1d103..48d2105 100644 (file)
              read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
              print *,msc(i,5),restok(i,5)
             enddo
-            ip(5)=0.2
+           ! ip(5)=0.2
            ! mp(5)=0.2
              pstok(5)=3.0
 !DIR$ NOUNROLL 
index c060c01..8a4020f 100644 (file)
@@ -24,7 +24,8 @@
 !
 !  ***  assess candidate step (***sol version 2.3)  ***
 !
-      integer :: liv, l,lv
+      integer :: liv, l
+      integer(kind=8) ::lv
       integer :: iv(liv)
       real(kind=8) :: v(lv)
 !
 !  ***  alg = 1 means regression constants.
 !  ***  alg = 2 means general unconstrained optimization constants.
 !
-      integer :: liv, l,lv
-      integer :: alg, iv(liv)
+      integer :: liv, l
+      integer(kind=8)::lv
+      integer :: iv(liv)
+      integer :: alg 
       real(kind=8) :: v(lv)
 !
 !el      external imdcon, vdflt
 !
 !  ***  parameter declarations  ***
 !
-      integer :: liv, lv, p
+      integer :: liv, p
+      integer(kind=8)::lv
       integer :: iv(liv)
       real(kind=8) :: d(p), g(p), v(lv), x(p)
 !
 !
 !  ***  alg = 1 for regression, alg = 2 for general unconstrained opt.
 !
-      integer :: alg, liv, lv, n
+      integer :: alg, liv, n
+      integer(kind=8) :: lv
       integer :: iv(liv)
       real(kind=8) :: d(n), v(lv)
 !
 !  ***  alg = 1 means regression constants.
 !  ***  alg = 2 means general unconstrained optimization constants.
 !
-      integer :: alg, l,lv
+      integer :: alg, l
+      integer(kind=8)::lv
       real(kind=8) :: v(lv)
 !/+
 !el      real(kind=8) :: dmax1
 !  ***  minimize general unconstrained objective function using   ***
 !  ***  (analytic) gradient and hessian provided by the caller.   ***
 !
-      integer :: liv, lv, n
-      integer :: iv(liv), uiparm(1)
+      integer :: liv, n
+      integer(kind=8):: lv
+      integer:: iv(liv)
+      integer :: uiparm(1)
       real(kind=8) :: d(n), x(n), v(lv), urparm(1)
       real(kind=8),external :: ufparm
 !     dimension v(78 + n*(n+12)), uiparm(*), urparm(*)
 
 !  ***  parameter declarations  ***
 !
-      integer :: lh, liv, lv, n
-      integer :: iv(liv)
+      integer :: lh, liv, n
+      integer(kind=8)::lv
+      integer ::  iv(liv)
       real(kind=8) :: d(n), fx, g(n), h(lh), v(lv), x(n)
 !
 !--------------------------  parameter usage  --------------------------
 !
 !  ***  parameter declarations  ***
 !
-      integer :: liv, lv, n
+      integer :: liv, n
+      integer(kind=8) ::lv
       integer :: iv(liv)
       real(kind=8) :: d(n), hdiag(n), v(lv)
 !
       logical :: not_done,change,reduce 
 !el      common /przechowalnia/ v
 !el local variables
-      integer :: iretcode,nfun,lv,nvar_restr,idum(1),j
+      integer :: iretcode,nfun,nvar_restr,idum(1),j
+      integer (kind=8):: lv
       real(kind=8) :: etot,rdum(1)     !,fdum
 
       lv=(77+(6*nres)*(6*nres+17)/2)   !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres)
 !      include 'COMMON.GEO'
 !      include 'COMMON.MINIM'
 !      include 'COMMON.CHAIN'
-      integer :: iretcode,nfun,k,i,j,lv,idum(1)
+      integer :: iretcode,nfun,k,i,j,idum(1)
+      integer(kind=8) :: lv
       integer,dimension(liv) :: iv                                               
       real(kind=8) :: minval   !,v(1:77+(6*nres)*(6*nres+17)/2)        !(1:lv)
       real(kind=8),dimension(6*nres) :: x,d,xx !(maxvar) (maxvar=6*maxres)
       real(kind=8) :: g(6*nres),f1,etot,rdum(1)        !,fdum
 
       lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres)
-
+      print *,"lv value",lv
       if (.not. allocated(v)) allocate(v(1:lv))
 
       call deflt(2,iv,liv,lv,v)                                         
+      print *,"after delft"
 ! 12 means fresh start, dont call deflt                                 
       iv(1)=12                                                          
 ! max num of fun calls                                                  
 !el      real(kind=8) :: v(1:77+(6*nres)*(6*nres+17)/2+1)                    
       integer,dimension(6) :: indx
       integer,dimension(liv) :: iv                                               
-      integer :: lv,idum(1),nf !
+      integer :: idum(1),nf    !
+      integer(kind=8) :: lv
       real(kind=8) :: rdum(1)
       real(kind=8) :: przes(3),obrot(3,3),eee
       logical :: non_conv
 !      external func_restr1,grad_restr1
       logical :: not_done,change,reduce 
 !el      common /przechowalnia/ v
-
-      integer :: nvar_restr,lv,i,j
+      integer(kind=8):: lv
+      integer :: nvar_restr,i,j
       integer :: idum(1)
       real(kind=8) :: rdum(1),etot     !,fdum
 
 !  ***  analytic gradient and hessian approx. from secant update  ***
 !
 !      use control
-      integer :: n, liv, lv
-      integer :: iv(liv), uiparm(1)
+      integer :: n, liv
+      integer(kind=8)::lv
+      integer:: iv(liv)
+      integer ::  uiparm(1)
       real(kind=8) :: d(n), x(n), v(lv), urparm(1)
       real(kind=8),external :: ufparm !funtion name as an argument
 
 ! sumit... reverse-communication routine that carries out sumsl algo-
 !             rithm.
 !
-      integer :: g1, iv1, nf
+      integer ::  iv1, nf
       real(kind=8) :: f
+      integer(kind=8)::g1
 !
 !  ***  subscripts for iv   ***
 !
       if (iv1 .eq. 12) iv(1) = 13
       go to 20
 !
- 10   g1 = iv(g)
+ 10   g1 = iv(g) 
+      print *,"my new g1",g1
 !elwrite(iout,*) "in sumsl go to 10"
 
 !
       subroutine sumit(d,fx,g,iv,liv,lv,n,v,x)
       
       use control, only:stopx
+      use MD_data,  only: itime_mat
 !
 !  ***  carry out sumsl (unconstrained minimization) iterations, using
 !  ***  double-dogleg/bfgs steps.
 !
 !  ***  parameter declarations  ***
 !
-      integer :: liv, lv, n
+      integer :: liv, n
+      integer(kind=8) :: lv
       integer :: iv(liv)
       real(kind=8) :: d(n), fx, g(n), v(lv), x(n)
 !
 !
  80   call itsum(d, g, iv, liv, lv, n, v, x)
  90   k = iv(niter)
+      itime_mat=k
+!      imat_update=k
       if (k .lt. iv(mxiter)) go to 100
          iv(1) = 10
          go to 300
 !
 !  ***  parameter declarations  ***
 !
-      integer :: lv, n
+      integer ::  n
+      integer(kind=8) :: lv
       real(kind=8) :: dig(n), nwtstp(n), step(n), v(lv)
 !
 !  ***  purpose  ***
index c53e992..e6c19b4 100644 (file)
 #endif
       print *,'Start MD'
       call alloc_MD_arrays
-      call alloc_compare_arrays
       print *,'After MD alloc'
+      call alloc_compare_arrays
+      print *,'After compare alloc'
       if (me.eq.king .or. .not. out1file) &
          write (iout,*) "Calling chainbuild"
       if (extconf) then