Adam's update
authorCezary Czaplewski <czarek@chem.univ.gda.pl>
Thu, 26 Mar 2020 17:49:14 +0000 (18:49 +0100)
committerCezary Czaplewski <czarek@chem.univ.gda.pl>
Thu, 26 Mar 2020 17:49:14 +0000 (18:49 +0100)
source/unres/src-HCD-5D/COMMON.INTERACT
source/unres/src-HCD-5D/COMMON.MD
source/unres/src-HCD-5D/MD_A-MTS.F
source/unres/src-HCD-5D/Makefile_MPICH_ifort-okeanos
source/unres/src-HCD-5D/energy_p_new_barrier.F
source/unres/src-HCD-5D/energy_p_new_barrier.F.safe
source/unres/src-HCD-5D/initialize_p.F
source/unres/src-HCD-5D/lagrangian_lesyng.F
source/unres/src-HCD-5D/make_xx_list.F [new file with mode: 0644]
source/unres/src-HCD-5D/readrtns_CSA.F
source/unres/src-HCD-5D/stochfric.F

index 14d92ef..3440239 100644 (file)
      & iscpstart(maxres,maxint_gr),iscpend(maxres,maxint_gr),
      & iatsc_s,iatsc_e,iatel_s,iatel_e,iatel_s_vdw,iatel_e_vdw,
      & iatscp_s,iatscp_e,ispp,iscp
+C 3/26/20 Interaction lists
+      integer newcontlisti(200*maxres),newcontlistj(200*maxres),
+     & newcontlistppi(200*maxres),newcontlistppj(200*maxres),
+     & newcontlistscpi(200*maxres),newcontlistscpj(200*maxres),
+     & g_listscsc_start,g_listscsc_end,g_listpp_start,g_listpp_end,
+     & g_listscp_start,g_listscp_end
+      common /interact_list/newcontlisti,newcontlistj,g_listscsc_start,
+     & g_listscsc_end,newcontlistppi,newcontlistppj,g_listpp_start,
+     & g_listpp_end,newcontlistscpi,newcontlistscpj,g_listscp_start,
+     & g_listscp_end
 C 12/1/95 Array EPS included in the COMMON block.
       double precision eps,epslip,sigma,sigmaii,rs0,chi,chip,alp,
      & sigma0,sigii,
index 6988bd8..cea18eb 100644 (file)
@@ -3,12 +3,12 @@
      & dvmax,damax,edriftmax
       integer n_timestep,ntwx,ntwe,lang,count_reset_moment,
      & count_reset_vel,ntime_split,ntime_split0,
-     & maxtime_split
+     & maxtime_split,itime_mat,imatupdate
       logical large,print_compon,tbf,rest,reset_moment,reset_vel,
      & rattle,mdpdb,RESPA,preminim
       common /mdpar/ v_ini,d_time,d_time0,t_bath,
      & tau_bath,dvmax,damax,n_timestep,mdpdb,
-     & ntime_split,ntime_split0,maxtime_split,
+     & ntime_split,ntime_split0,maxtime_split,itime_mat,imatupdate,
      & ntwx,ntwe,lang,large,print_compon,tbf,rest,preminim,
      & reset_moment,reset_vel,count_reset_moment,count_reset_vel,
      & rattle,RESPA
index 08852ba..4498b4a 100644 (file)
@@ -166,14 +166,14 @@ c   Entering the MD loop
      &      .and. mod(itime,count_reset_vel).eq.0) then
           call random_vel
           write(iout,'(a,f20.2)') 
-     &     "Velocities reset to random values, time",totT      
+     &     "Velocities reset to random values, time",totT
           do i=0,2*nres
             do j=1,3
               d_t_old(j,i)=d_t(j,i)
             enddo
           enddo
         endif
-               if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
+        if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
           call inertia_tensor  
           call vcm_vel(vcm)
           do j=1,3
@@ -182,7 +182,7 @@ c   Entering the MD loop
           call kinetic(EK)
           kinetic_T=2.0d0/(dimen3*Rb)*EK
           scalfac=dsqrt(T_bath/kinetic_T)
-          write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT      
+          write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
           do i=0,2*nres
             do j=1,3
               d_t_old(j,i)=scalfac*d_t(j,i)
@@ -209,6 +209,7 @@ c Variable time step algorithm.
           stop
 #endif
         endif
+        itime_mat=itime
         if (ntwe.ne.0) then
          if (mod(itime,ntwe).eq.0) then
            call statout(itime)
@@ -1811,10 +1812,11 @@ c      rest2name = prefix(:ilen(prefix))//'.rst'
          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
      &   (d_t(j,i+nres),j=1,3)
        enddo
+      endif
+      if (lang.eq.0) then
 c  Zeroing the total angular momentum of the system
        write(iout,*) "Calling the zero-angular 
      & momentum subroutine"
-      endif
       call inertia_tensor  
 c  Getting the potential energy and forces and velocities and accelerations
       call vcm_vel(vcm)
@@ -1829,8 +1831,14 @@ c Removing the velocity of the center of mass
       if(me.eq.king.or..not.out1file)then
        write (iout,*) "vcm right after adjustment:"
        write (iout,*) (vcm(j),j=1,3) 
+       write (iout,*) "Initial velocities after adjustment"
+       do i=0,nres
+         write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
+     &   (d_t(j,i+nres),j=1,3)
+       enddo
+       call flush(iout)
+      endif
       endif
-      call flush(iout)
       write (iout,*) "init_MD before initial structure REST ",rest
       if (.not.rest) then
   122   continue
@@ -2050,6 +2058,7 @@ c       write(iout,*) (potEcomp(i),i=0,n_ene)
 c      write (iout,*) "PotE-homology",potE
       totE=EK+potE
       itime=0
+      itime_mat=itime
       if (ntwe.ne.0) call statout(itime)
       if(me.eq.king.or..not.out1file)
      &  write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
@@ -2351,6 +2360,9 @@ c        write (iout,*) "ichain",ichain," innt",innt," inct",inct
 !
 #define WLOS
 #ifdef WLOS
+      if (nnt.eq.1) then
+        d_t(:,0)=d_t(:,1)
+      endif
       do i=1,nres
         if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then
           do j=1,3
@@ -2363,11 +2375,17 @@ c        write (iout,*) "ichain",ichain," innt",innt," inct",inct
           enddo
         end if
       enddo
+      d_t(:,nres)=0.0d0
       d_t(:,nct)=0.0d0
+      d_t(:,2*nres)=0.0d0
+      if (nnt.gt.1) then
+        d_t(:,0)=d_t(:,1)
+        d_t(:,1)=0.0d0
+      endif
 c      d_a(:,0)=d_a(:,1)
 c      d_a(:,1)=0.0d0
 c      write (iout,*) "Shifting accelerations"
-      do ichain=1,nchain
+      do ichain=2,nchain
 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))
index 5485424..32a0dec 100644 (file)
@@ -38,7 +38,7 @@ object = unres.o arcos.o cartprint.o chainbuild.o convert.o initialize_p.o \
         matmult.o readrtns_CSA.o parmread.o gen_rand_conf.o printmat.o map.o \
         pinorm.o randgens.o rescode.o intcor.o timing.o misc.o \
         cart2intgrad.o checkder_p.o contact_cp econstr_local.o econstr_qlike.o \
-       econstrq-PMF.o PMFprocess.o energy_p_new_barrier.o \
+       econstrq-PMF.o PMFprocess.o energy_p_new_barrier.o make_xx_list \
        energy_p_new-sep_barrier.o gradient_p.o minimize_p.o sumsld.o \
         cored.o rmdd.o geomout.o readpdb.o regularize.o thread.o fitsq.o mcm.o \
         mc.o bond_move.o refsys.o check_sc_distr.o check_bond.o contact.o \
@@ -95,7 +95,7 @@ NEWCORR: ${object} xdrf/libxdrf.a
 
 NEWCORR5D: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DAMD64 -DUNRES -DISNAN -DMP -DMPI \
        -DSPLITELE -DLANG0 -DNEWCORR -DCORRCD -DFIVEDIAG -DLBFGS #-DMYGAUSS #-DTIMING
-NEWCORR5D: BIN = ~/bin/unres_ifort_MPICH-okeanos_SC-HCD5.exe
+NEWCORR5D: BIN = ~/bin/unres_ifort_MPICH-okeanos_SC-HCD5-l.exe
 NEWCORR5D: ${object_lbfgs} ${object} fdisy.o fdiag.o machpd.o kinetic_CASC.o xdrf/libxdrf.a
        gcc -o compinfo compinfo.c
        ./compinfo | true
index 2a588bd..ef19809 100644 (file)
@@ -30,6 +30,7 @@ c      include 'COMMON.MD'
       include 'COMMON.SPLITELE'
       include 'COMMON.TORCNSTR'
       include 'COMMON.SAXS'
+      include 'COMMON.MD'
       double precision evdw,evdw1,evdw2,evdw2_14,ees,eel_loc,
      & eello_turn3,eello_turn4,edfadis,estr,ehpb,ebe,ethetacnstr,
      & escloc,etors,edihcnstr,etors_d,esccor,ecorr,ecorr5,ecorr6,eturn6,
@@ -117,6 +118,12 @@ c        call chainbuild_cart
       edfanei=0.0d0
       edfabet=0.0d0
 #endif
+      if (nfgtasks.gt.1) then
+        call MPI_Bcast(itime_mat,1,MPI_INT,king,FG_COMM,IERROR)
+      endif
+      if (mod(itime_mat,imatupdate).eq.0) call make_SCp_inter_list
+      if (mod(itime_mat,imatupdate).eq.0) call make_SCSC_inter_list
+      if (mod(itime_mat,imatupdate).eq.0) call make_pp_inter_list
 c      print *,'Processor',myrank,' calling etotal ipot=',ipot
 c      print *,'Processor',myrank,' nnt=',nnt,' nct=',nct
 #else
 #endif
       double precision gg(3)
       double precision evdw,evdwij
-      integer i,j,k,itypi,itypj,itypi1,num_conti,iint
+      integer i,j,k,itypi,itypj,itypi1,num_conti,iint,icont
       double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
      & sigij,r0ij,rcut,sqrij,sss1,sssgrad1
       double precision fcont,fprimcont
       double precision sscale,sscagrad
 c      write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
-      do i=iatsc_s,iatsc_e
+c      do i=iatsc_s,iatsc_e
+      do icont=g_listscsc_start,g_listscsc_end
+        i=newcontlisti(icont)
+        j=newcontlistj(icont)
         itypi=iabs(itype(i))
         if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
@@ -1479,10 +1489,10 @@ C Change 12/1/95
 C
 C Calculate SC interaction energy.
 C
-        do iint=1,nint_gr(i)
+c        do iint=1,nint_gr(i)
 cd        write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
 cd   &                  'iend=',iend(i,iint)
-          do j=istart(i,iint),iend(i,iint)
+c          do j=istart(i,iint),iend(i,iint)
             itypj=iabs(itype(j)) 
             if (itypj.eq.ntyp1) cycle
             xj=c(1,nres+j)-xi
@@ -1587,8 +1597,8 @@ cd   &           i,j,(gacont(kk,num_conti,i),kk=1,3)
               endif
             endif
 #endif
-          enddo      ! j
-        enddo        ! iint
+c          enddo      ! j
+c        enddo        ! iint
 C Change 12/1/95
 #ifdef FOURBODY
         num_cont(i)=num_conti
       include 'COMMON.SPLITELE'
       double precision gg(3)
       double precision evdw,evdwij
-      integer i,j,k,itypi,itypj,itypi1,iint
+      integer i,j,k,itypi,itypj,itypi1,iint,icont
       double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
      & fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1
       logical scheck
       double precision sscale,sscagrad
 c     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
-      do i=iatsc_s,iatsc_e
+c      do i=iatsc_s,iatsc_e
+      do icont=g_listscsc_start,g_listscsc_end
+        i=newcontlisti(icont)
+        j=newcontlistj(icont)
         itypi=iabs(itype(i))
         if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
@@ -1647,8 +1660,8 @@ c     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
 C
 C Calculate SC interaction energy.
 C
-        do iint=1,nint_gr(i)
-          do j=istart(i,iint),iend(i,iint)
+c        do iint=1,nint_gr(i)
+c          do j=istart(i,iint),iend(i,iint)
             itypj=iabs(itype(j))
             if (itypj.eq.ntyp1) cycle
             xj=c(1,nres+j)-xi
@@ -1695,8 +1708,8 @@ cgrad              do l=1,3
 cgrad                gvdwc(l,k)=gvdwc(l,k)+gg(l)
 cgrad              enddo
 cgrad            enddo
-          enddo      ! j
-        enddo        ! iint
+c          enddo      ! j
+c        enddo        ! iint
       enddo          ! i
       do i=1,nct
         do j=1,3
@@ -1727,7 +1740,7 @@ C
       integer icall
       common /srutu/ icall
       double precision evdw
-      integer itypi,itypj,itypi1,iint,ind
+      integer itypi,itypj,itypi1,iint,ind,icont
       double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi,
      & sss1,sssgrad1
       double precision sscale,sscagrad
@@ -1742,7 +1755,10 @@ c     else
         lprn=.false.
 c     endif
       ind=0
-      do i=iatsc_s,iatsc_e
+c      do i=iatsc_s,iatsc_e 
+      do icont=g_listscsc_start,g_listscsc_end
+        i=newcontlisti(icont)
+        j=newcontlistj(icont)
         itypi=iabs(itype(i))
         if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
@@ -1757,8 +1773,8 @@ c        dsci_inv=dsc_inv(itypi)
 C
 C Calculate SC interaction energy.
 C
-        do iint=1,nint_gr(i)
-          do j=istart(i,iint),iend(i,iint)
+c        do iint=1,nint_gr(i)
+c          do j=istart(i,iint),iend(i,iint)
             ind=ind+1
             itypj=iabs(itype(j))
             if (itypj.eq.ntyp1) cycle
@@ -1835,8 +1851,8 @@ C Calculate radial part of the gradient
 C Calculate the angular part of the gradient and sum add the contributions
 C to the appropriate components of the Cartesian gradient.
             call sc_grad
-          enddo      ! j
-        enddo        ! iint
+!          enddo      ! j
+!        enddo        ! iint
       enddo          ! i
 c     stop
       return
@@ -1864,7 +1880,7 @@ C
       logical lprn
       integer xshift,yshift,zshift,subchap
       double precision evdw
-      integer itypi,itypj,itypi1,iint,ind
+      integer itypi,itypj,itypi1,iint,ind,icont
       double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
       double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
      & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
@@ -1882,7 +1898,10 @@ C we have the original box)
 C      do xshift=-1,1
 C      do yshift=-1,1
 C      do zshift=-1,1
-      do i=iatsc_s,iatsc_e
+c      do i=iatsc_s,iatsc_e
+      do icont=g_listscsc_start,g_listscsc_end
+        i=newcontlisti(icont)
+        j=newcontlistj(icont)
         itypi=iabs(itype(i))
         if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
@@ -1963,8 +1982,8 @@ c        write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
 C
 C Calculate SC interaction energy.
 C
-        do iint=1,nint_gr(i)
-          do j=istart(i,iint),iend(i,iint)
+c        do iint=1,nint_gr(i)
+c          do j=istart(i,iint),iend(i,iint)
             IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
 
 c              write(iout,*) "PRZED ZWYKLE", evdwij
@@ -2211,8 +2230,8 @@ C Calculate angular part of the gradient.
 c            call sc_grad_scale(sss)
             call sc_grad
             ENDIF    ! dyn_ss            
-          enddo      ! j
-        enddo        ! iint
+c          enddo      ! j
+c        enddo        ! iint
       enddo          ! i
 C      enddo          ! zshift
 C      enddo          ! yshift
@@ -2244,7 +2263,7 @@ C
       common /srutu/ icall
       logical lprn
       double precision evdw
-      integer itypi,itypj,itypi1,iint,ind
+      integer itypi,itypj,itypi1,iint,ind,icont
       double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,r0ij,
      & xi,yi,zi,fac_augm,e_augm
       double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
@@ -2257,7 +2276,10 @@ c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
       lprn=.false.
 c     if (icall.eq.0) lprn=.true.
       ind=0
-      do i=iatsc_s,iatsc_e
+c      do i=iatsc_s,iatsc_e
+      do icont=g_listscsc_start,g_listscsc_end
+        i=newcontlisti(icont)
+        j=newcontlistj(icont)
         itypi=iabs(itype(i))
         if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
@@ -2307,8 +2329,8 @@ c        dsci_inv=dsc_inv(itypi)
 C
 C Calculate SC interaction energy.
 C
-        do iint=1,nint_gr(i)
-          do j=istart(i,iint),iend(i,iint)
+c        do iint=1,nint_gr(i)
+c          do j=istart(i,iint),iend(i,iint)
             ind=ind+1
             itypj=iabs(itype(j))
             if (itypj.eq.ntyp1) cycle
@@ -2460,8 +2482,8 @@ C Calculate the radial part of the gradient
 C Calculate angular part of the gradient.
 c            call sc_grad_scale(sss)
             call sc_grad
-          enddo      ! j
-        enddo        ! iint
+c          enddo      ! j
+c        enddo        ! iint
       enddo          ! i
       end
 C-----------------------------------------------------------------------------
@@ -2612,7 +2634,10 @@ c      include 'COMMON.CONTACTS'
       dimension gg(3)
 cd    print *,'Entering Esoft_sphere nnt=',nnt,' nct=',nct
       evdw=0.0D0
-      do i=iatsc_s,iatsc_e
+c      do i=iatsc_s,iatsc_e
+      do icont=g_listscsc_start,g_listscsc_end
+        i=newcontlisti(icont)
+        j=newcontlistj(icont)
         itypi=iabs(itype(i))
         if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
@@ -2622,10 +2647,10 @@ cd    print *,'Entering Esoft_sphere nnt=',nnt,' nct=',nct
 C
 C Calculate SC interaction energy.
 C
-        do iint=1,nint_gr(i)
+c        do iint=1,nint_gr(i)
 cd        write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
 cd   &                  'iend=',iend(i,iint)
-          do j=istart(i,iint),iend(i,iint)
+c          do j=istart(i,iint),iend(i,iint)
             itypj=iabs(itype(j))
             if (itypj.eq.ntyp1) cycle
             xj=c(1,nres+j)-xi
@@ -2661,8 +2686,8 @@ cgrad              do l=1,3
 cgrad                gvdwc(l,k)=gvdwc(l,k)+gg(l)
 cgrad              enddo
 cgrad            enddo
-          enddo ! j
-        enddo ! iint
+c          enddo ! j
+c        enddo ! iint
       enddo ! i
       return
       end
@@ -3843,7 +3868,10 @@ c
 c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
 c
 CTU KURWA
-      do i=iatel_s,iatel_e
+c      do i=iatel_s,iatel_e
+       do icont=g_listpp_start,g_listpp_end
+        i=newcontlistppi(icont)
+        j=newcontlistppj(icont)
 C        do i=75,75
 c        if (i.le.1) cycle
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
@@ -3904,7 +3932,7 @@ c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
         num_conti=num_cont_hb(i)
 #endif
 C I TU KURWA
-        do j=ielstart(i),ielend(i)
+c        do j=ielstart(i),ielend(i)
 C          do j=16,17
 C          write (iout,*) i,j
 C         if (j.le.1) cycle
@@ -3917,7 +3945,7 @@ c     & .or.itype(j+2).eq.ntyp1
 c     & .or.itype(j-1).eq.ntyp1
      &) cycle
           call eelecij(i,j,ees,evdw1,eel_loc)
-        enddo ! j
+c        enddo ! j
 #ifdef FOURBODY
         num_cont_hb(i)=num_conti
 #endif
@@ -5546,7 +5574,10 @@ cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
 C      do xshift=-1,1
 C      do yshift=-1,1
 C      do zshift=-1,1
-      do i=iatscp_s,iatscp_e
+c      do i=iatscp_s,iatscp_e
+      do icont=g_listscp_start,g_listscp_end
+        i=newcontlistscpi(icont)
+        j=newcontlistscpj(icont)
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
         iteli=itel(i)
         xi=0.5D0*(c(1,i)+c(1,i+1))
@@ -5586,9 +5617,9 @@ c        endif
 C          xi=xi+xshift*boxxsize
 C          yi=yi+yshift*boxysize
 C          zi=zi+zshift*boxzsize
-        do iint=1,nscp_gr(i)
+c        do iint=1,nscp_gr(i)
 
-        do j=iscpstart(i,iint),iscpend(i,iint)
+c        do j=iscpstart(i,iint),iscpend(i,iint)
           if (itype(j).eq.ntyp1) cycle
           itypj=iabs(itype(j))
 C Uncomment following three lines for SC-p interactions
@@ -5711,9 +5742,9 @@ cgrad          enddo
             gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
             gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
           enddo
-        enddo
+c        enddo
 
-        enddo ! iint
+c        enddo ! iint
       enddo ! i
 C      enddo !zshift
 C      enddo !yshift
@@ -5741,7 +5772,7 @@ C
       include 'COMMON.SPLITELE'
       integer xshift,yshift,zshift
       double precision ggg(3)
-      integer i,iint,j,k,iteli,itypj,subchap
+      integer i,iint,j,k,iteli,itypj,subchap,icont
       double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1,
      & fac,e1,e2,rij
       double precision evdw2,evdw2_14,evdwij
@@ -5757,7 +5788,10 @@ C      do xshift=-1,1
 C      do yshift=-1,1
 C      do zshift=-1,1
       if (energy_dec) write (iout,*) "escp:",r_cut_int,rlamb
-      do i=iatscp_s,iatscp_e
+c      do i=iatscp_s,iatscp_e
+      do icont=g_listscp_start,g_listscp_end
+        i=newcontlistscpi(icont)
+        j=newcontlistscpj(icont)
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
         iteli=itel(i)
         xi=0.5D0*(c(1,i)+c(1,i+1))
@@ -5800,9 +5834,9 @@ c        if ((zi.gt.((zshift+0.5d0)*boxzsize)).or.
 c     &       (zi.lt.((zshift-0.5d0)*boxzsize))) then
 c        go to 136
 c        endif
-        do iint=1,nscp_gr(i)
+c        do iint=1,nscp_gr(i)
 
-        do j=iscpstart(i,iint),iscpend(i,iint)
+c        do j=iscpstart(i,iint),iscpend(i,iint)
           itypj=iabs(itype(j))
           if (itypj.eq.ntyp1) cycle
 C Uncomment following three lines for SC-p interactions
@@ -5937,9 +5971,9 @@ cgrad          enddo
             gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
           enddo
 c        endif !endif for sscale cutoff
-        enddo ! j
+c        enddo ! j
 
-        enddo ! iint
+c        enddo ! iint
       enddo ! i
 c      enddo !zshift
 c      enddo !yshift
index ae8e449..2a588bd 100644 (file)
@@ -66,7 +66,8 @@ C FG slaves as WEIGHTS array.
           weights_(17)=wbond
           weights_(18)=scal14
           weights_(21)=wsccor
-          weights_(22)=wtube
+          weights_(22)=wliptran
+          weights_(25)=wtube
           weights_(26)=wsaxs
           weights_(28)=wdfa_dist
           weights_(29)=wdfa_tor
@@ -98,7 +99,8 @@ C FG slaves receive the WEIGHTS array
           wbond=weights(17)
           scal14=weights(18)
           wsccor=weights(21)
-          wtube=weights(22)
+          wliptran=weights(22)
+          wtube=weights(25)
           wsaxs=weights(26)
           wdfa_dist=weights_(28)
           wdfa_tor=weights_(29)
@@ -387,6 +389,8 @@ C based on partition function
 C      print *,"przed lipidami"
       if (wliptran.gt.0) then
         call Eliptransfer(eliptran)
+      else
+        eliptran=0.0d0
       endif
 C      print *,"za lipidami"
       if (AFMlog.gt.0) then
@@ -1449,6 +1453,7 @@ C
       include 'COMMON.SBRIDGE'
       include 'COMMON.NAMES'
       include 'COMMON.IOUNITS'
+      include 'COMMON.SPLITELE'
 #ifdef FOURBODY
       include 'COMMON.CONTACTS'
       include 'COMMON.CONTMAT'
@@ -1457,8 +1462,9 @@ C
       double precision evdw,evdwij
       integer i,j,k,itypi,itypj,itypi1,num_conti,iint
       double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
-     & sigij,r0ij,rcut
+     & sigij,r0ij,rcut,sqrij,sss1,sssgrad1
       double precision fcont,fprimcont
+      double precision sscale,sscagrad
 c      write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
       do i=iatsc_s,iatsc_e
@@ -1485,6 +1491,11 @@ cd   &                  'iend=',iend(i,iint)
 C Change 12/1/95 to calculate four-body interactions
             rij=xj*xj+yj*yj+zj*zj
             rrij=1.0D0/rij
+            sqrij=dsqrt(rij)
+            sss1=sscale(sqrij,r_cut_int)
+            if (sss1.eq.0.0d0) cycle
+            sssgrad1=sscagrad(sqrij,r_cut_int)
+            
 c           write (iout,*)'i=',i,' j=',j,' itypi=',itypi,' itypj=',itypj
             eps0ij=eps(itypi,itypj)
             fac=rrij**expon2
@@ -1498,11 +1509,12 @@ cd          write (iout,'(2(a3,i3,2x),6(1pd12.4)/2(3(1pd12.4),5x)/)')
 cd   &        restyp(itypi),i,restyp(itypj),j,a(itypi,itypj),
 cd   &        bb(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,epsi,sigm,
 cd   &        (c(k,i),k=1,3),(c(k,j),k=1,3)
-            evdw=evdw+evdwij
+            evdw=evdw+sss1*evdwij
 C 
 C Calculate the components of the gradient in DC and X
 C
-            fac=-rrij*(e1+evdwij)
+            fac=-rrij*(e1+evdwij)*sss1
+     &          +evdwij*sssgrad1/sqrij/expon
             gg(1)=xj*fac
             gg(2)=yj*fac
             gg(3)=zj*fac
       include 'COMMON.INTERACT'
       include 'COMMON.IOUNITS'
       include 'COMMON.NAMES'
+      include 'COMMON.SPLITELE'
       double precision gg(3)
       double precision evdw,evdwij
       integer i,j,k,itypi,itypj,itypi1,iint
       double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
-     & fac_augm,e_augm,r_inv_ij,r_shift_inv
+     & fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1
       logical scheck
+      double precision sscale,sscagrad
 c     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
       do i=iatsc_s,iatsc_e
@@ -1645,6 +1659,9 @@ C
             e_augm=augm(itypi,itypj)*fac_augm
             r_inv_ij=dsqrt(rrij)
             rij=1.0D0/r_inv_ij 
+            sss1=sscale(rij,r_cut_int)
+            if (sss1.eq.0.0d0) cycle
+            sssgrad1=sscagrad(rij,r_cut_int)
             r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
             fac=r_shift_inv**expon
 C have you changed here?
@@ -1658,11 +1675,12 @@ cd   &        restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
 cd   &        bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
 cd   &        sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
 cd   &        (c(k,i),k=1,3),(c(k,j),k=1,3)
-            evdw=evdw+evdwij
+            evdw=evdw+evdwij*sss1
 C 
 C Calculate the components of the gradient in DC and X
 C
             fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
+     &          +evdwij*sssgrad1*r_inv_ij/expon
             gg(1)=xj*fac
             gg(2)=yj*fac
             gg(3)=zj*fac
       include 'COMMON.INTERACT'
       include 'COMMON.IOUNITS'
       include 'COMMON.CALC'
+      include 'COMMON.SPLITELE'
       integer icall
       common /srutu/ icall
       double precision evdw
       integer itypi,itypj,itypi1,iint,ind
-      double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
+      double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi,
+     & sss1,sssgrad1
+      double precision sscale,sscagrad
 c     double precision rrsave(maxdim)
       logical lprn
       evdw=0.0D0
@@ -1775,6 +1796,9 @@ cd          else
 cd            rrij=rrsave(ind)
 cd          endif
             rij=dsqrt(rrij)
+            sss1=sscale(1.0d0/rij,r_cut_int)
+            if (sss1.eq.0.0d0) cycle
+            sssgrad1=sscagrad(1.0d0/rij,r_cut_int)
 C Calculate the angle-dependent terms of energy & contributions to derivatives.
             call sc_angular
 C Calculate whole angle-dependent part of epsilon and contributions
@@ -1787,7 +1811,7 @@ C have you changed here?
             eps2der=evdwij*eps3rt
             eps3der=evdwij*eps2rt
             evdwij=evdwij*eps2rt*eps3rt
-            evdw=evdw+evdwij
+            evdw=evdw+sss1*evdwij
             if (lprn) then
             sigm=dabs(aa/bb)**(1.0D0/6.0D0)
             epsi=bb**2/aa
@@ -1803,6 +1827,7 @@ C Calculate gradient components.
             fac=-expon*(e1+evdwij)
             sigder=fac/sigsq
             fac=rrij*fac
+     &          +evdwij*sssgrad1/sss1*rij
 C Calculate radial part of the gradient
             gg(1)=xj*fac
             gg(2)=yj*fac
@@ -2108,11 +2133,10 @@ c     &       dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
             rij=dsqrt(rrij)
             sss=sscale(1.0d0/rij,r_cut_int)
-            sssgrad=sscagrad(1.0d0/rij,r_cut_int)
-             
 c            write (iout,'(a7,4f8.3)') 
 c    &      "ssscale",sss,((1.0d0/rij)/sigma(itypi,itypj)),r_cut,rlamb
-            if (sss.gt.0.0d0) then
+            if (sss.eq.0.0d0) cycle
+            sssgrad=sscagrad(1.0d0/rij,r_cut_int)
 C Calculate angle-dependent terms of energy and contributions to their
 C derivatives.
             call sc_angular
@@ -2159,8 +2183,8 @@ c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
      &        evdwij
             endif
 
-            if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') 
-     &                        'evdw',i,j,evdwij
+            if (energy_dec) write (iout,'(a,2i5,3f10.5)') 
+     &                    'r sss evdw',i,j,rij,sss,evdwij
 
 C Calculate gradient components.
             e1=e1*eps1*eps2rt**2*eps3rt**2
@@ -2169,13 +2193,13 @@ C Calculate gradient components.
             fac=rij*fac
 c            print '(2i4,6f8.4)',i,j,sss,sssgrad*
 c     &      evdwij,fac,sigma(itypi,itypj),expon
-            fac=fac+evdwij/sss*sssgrad*rij
+            fac=fac+evdwij*sssgrad/sss*rij
 c            fac=0.0d0
 C Calculate the radial part of the gradient
             gg_lipi(3)=eps1*(eps2rt*eps2rt)
-     &*(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip*
-     & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
-     &+faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj)))
+     &       *(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip*
+     &        (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
+     &       +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj)))
             gg_lipj(3)=ssgradlipj*gg_lipi(3)
             gg_lipi(3)=gg_lipi(3)*ssgradlipi
 C            gg_lipi(3)=0.0d0
@@ -2184,8 +2208,8 @@ C            gg_lipj(3)=0.0d0
             gg(2)=yj*fac
             gg(3)=zj*fac
 C Calculate angular part of the gradient.
+c            call sc_grad_scale(sss)
             call sc_grad
-            endif
             ENDIF    ! dyn_ss            
           enddo      ! j
         enddo        ! iint
@@ -2214,6 +2238,7 @@ C
       include 'COMMON.INTERACT'
       include 'COMMON.IOUNITS'
       include 'COMMON.CALC'
+      include 'COMMON.SPLITELE'
       integer xshift,yshift,zshift,subchap
       integer icall
       common /srutu/ icall
@@ -2224,7 +2249,7 @@ C
      & xi,yi,zi,fac_augm,e_augm
       double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
      & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
-     & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
+     & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip,sssgrad1
       double precision dist,sscale,sscagrad,sscagradlip,sscalelip
       evdw=0.0D0
 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
@@ -2384,6 +2409,9 @@ C      write(iout,*) "tu,", i,j,aa,bb,aa_lip(itypi,itypj),sslipi,sslipj
             dzj=dc_norm(3,nres+j)
             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
             rij=dsqrt(rrij)
+            sss=sscale(1.0d0/rij,r_cut_int)
+            if (sss.eq.0.0d0) cycle
+            sssgrad=sscagrad(1.0d0/rij,r_cut_int)
 C Calculate angle-dependent terms of energy and contributions to their
 C derivatives.
             call sc_angular
@@ -2424,12 +2452,13 @@ C Calculate gradient components.
             fac=-expon*(e1+evdwij)*rij_shift
             sigder=fac*sigder
             fac=rij*fac-2*expon*rrij*e_augm
-            fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij
+            fac=fac+(evdwij+e_augm)*sssgrad/sss*rij
 C Calculate the radial part of the gradient
             gg(1)=xj*fac
             gg(2)=yj*fac
             gg(3)=zj*fac
 C Calculate angular part of the gradient.
+c            call sc_grad_scale(sss)
             call sc_grad
           enddo      ! j
         enddo        ! iint
@@ -3909,7 +3938,7 @@ cd      print *,"Processor",fg_rank," t_eelecij",t_eelecij
       end
 C-------------------------------------------------------------------------------
       subroutine eelecij(i,j,ees,evdw1,eel_loc)
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
 #ifdef MPI
       include "mpif.h"
@@ -3933,14 +3962,33 @@ C-------------------------------------------------------------------------------
       include 'COMMON.TIME1'
       include 'COMMON.SPLITELE'
       include 'COMMON.SHIELD'
-      dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
+      double precision ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
      &          erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
       double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
      &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),gmuij1(4),gmuji1(4),
      &    gmuij2(4),gmuji2(4)
+      double precision dxi,dyi,dzi
+      double precision dx_normi,dy_normi,dz_normi,aux
+      integer j1,j2,lll,num_conti
       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
      &    dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
      &    num_conti,j1,j2
+      integer k,i,j,iteli,itelj,kkk,l,kkll,m,isubchap,ilist,iresshield
+      double precision ael6i,rrmij,rmij,r0ij,fcont,fprimcont,ees0tmp
+      double precision ees,evdw1,eel_loc,aaa,bbb,ael3i
+      double precision dxj,dyj,dzj,dx_normj,dy_normj,dz_normj,xj,yj,zj,
+     &  rij,r3ij,r6ij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4,
+     &  evdwij,el1,el2,eesij,ees0ij,facvdw,facel,fac1,ecosa,
+     &  ecosb,ecosg,ury,urz,vry,vrz,facr,a22der,a23der,a32der,
+     &  a33der,eel_loc_ij,cosa4,wij,cosbg1,cosbg2,ees0pij,
+     &  ees0pij1,ees0mij,ees0mij1,fac3p,ees0mijp,ees0pijp,
+     &  ecosa1,ecosb1,ecosg1,ecosa2,ecosb2,ecosg2,ecosap,ecosbp,
+     &  ecosgp,ecosam,ecosbm,ecosgm,ghalf,rlocshield
+      double precision a22,a23,a32,a33,geel_loc_ij,geel_loc_ji
+      double precision dist_init,xj_safe,yj_safe,zj_safe,
+     &  xj_temp,yj_temp,zj_temp,dist_temp,xmedi,ymedi,zmedi
+      double precision sscale,sscagrad,scalar
+
 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
 #ifdef MOMENT
       double precision scal_el /1.0d0/
@@ -4044,8 +4092,9 @@ C        yj=yj-ymedi
 C        zj=zj-zmedi
           rij=xj*xj+yj*yj+zj*zj
 
-            sss=sscale(sqrt(rij),r_cut_int)
-            sssgrad=sscagrad(sqrt(rij),r_cut_int)
+          sss=sscale(dsqrt(rij),r_cut_int)
+          if (sss.eq.0.0d0) return
+          sssgrad=sscagrad(dsqrt(rij),r_cut_int)
 c            if (sss.gt.0.0d0) then  
           rrmij=1.0D0/rij
           rij=dsqrt(rij)
@@ -4080,7 +4129,7 @@ C          fac_shield(j)=0.6
           fac_shield(i)=1.0
           fac_shield(j)=1.0
           eesij=(el1+el2)
-          ees=ees+eesij
+          ees=ees+eesij*sss
           endif
           evdw1=evdw1+evdwij*sss
 cd          write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
@@ -4089,11 +4138,10 @@ cd     &      1.0D0/dsqrt(rrmij),evdwij,eesij,
 cd     &      xmedi,ymedi,zmedi,xj,yj,zj
 
           if (energy_dec) then 
-              write (iout,'(a6,2i5,0pf7.3,2i5,3e11.3)') 
-     &'evdw1',i,j,evdwij
-     &,iteli,itelj,aaa,evdw1,sss
-              write (iout,'(a6,2i5,0pf7.3,2f8.3)') 'ees',i,j,eesij,
-     &fac_shield(i),fac_shield(j)
+            write (iout,'(a6,2i5,0pf7.3,2i5,e11.3,3f10.5)') 
+     &        'evdw1',i,j,evdwij,iteli,itelj,aaa,evdw1,sss,rij
+            write (iout,'(a6,2i5,0pf7.3,2f8.3)') 'ees',i,j,eesij,
+     &        fac_shield(i),fac_shield(j)
           endif
 
 C
@@ -4110,9 +4158,10 @@ C
 *
 * Radial derivatives. First process both termini of the fragment (i,j)
 *
-          ggg(1)=facel*xj
-          ggg(2)=facel*yj
-          ggg(3)=facel*zj
+          aux=facel*sss+rmij*sssgrad*eesij
+          ggg(1)=aux*xj
+          ggg(2)=aux*yj
+          ggg(3)=aux*zj
           if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
      &  (shield_mode.gt.0)) then
 C          print *,i,j     
@@ -4146,10 +4195,10 @@ C              endif
            iresshield=shield_list(ilist,j)
            do k=1,3
            rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j)
-     &     *2.0
+     &     *2.0*sss
            gshieldx(k,iresshield)=gshieldx(k,iresshield)+
      &              rlocshield
-     & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0
+     & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0*sss
            gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
 
 C     & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)
@@ -4172,13 +4221,13 @@ C              endif
 
           do k=1,3
             gshieldc(k,i)=gshieldc(k,i)+
-     &              grad_shield(k,i)*eesij/fac_shield(i)*2.0
+     &              grad_shield(k,i)*eesij/fac_shield(i)*2.0*sss
             gshieldc(k,j)=gshieldc(k,j)+
-     &              grad_shield(k,j)*eesij/fac_shield(j)*2.0
+     &              grad_shield(k,j)*eesij/fac_shield(j)*2.0*sss
             gshieldc(k,i-1)=gshieldc(k,i-1)+
-     &              grad_shield(k,i)*eesij/fac_shield(i)*2.0
+     &              grad_shield(k,i)*eesij/fac_shield(i)*2.0*sss
             gshieldc(k,j-1)=gshieldc(k,j-1)+
-     &              grad_shield(k,j)*eesij/fac_shield(j)*2.0
+     &              grad_shield(k,j)*eesij/fac_shield(j)*2.0*sss
 
            enddo
            endif
@@ -4209,15 +4258,10 @@ cgrad            do l=1,3
 cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
 cgrad            enddo
 cgrad          enddo
-          if (sss.gt.0.0) then
-          ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj
-          ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj
-          ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj
-          else
-          ggg(1)=0.0
-          ggg(2)=0.0
-          ggg(3)=0.0
-          endif
+          facvdw=facvdw+sssgrad*rmij*evdwij
+          ggg(1)=facvdw*xj
+          ggg(2)=facvdw*yj
+          ggg(3)=facvdw*zj
 c          do k=1,3
 c            ghalf=0.5D0*ggg(k)
 c            gvdwpp(k,i)=gvdwpp(k,i)+ghalf
@@ -4238,10 +4282,11 @@ cgrad            enddo
 cgrad          enddo
 #else
 C MARYSIA
-          facvdw=(ev1+evdwij)*sss
+          facvdw=(ev1+evdwij)
           facel=(el1+eesij)
           fac1=fac
-          fac=-3*rrmij*(facvdw+facvdw+facel)
+          fac=-3*rrmij*(facvdw+facvdw+facel)*sss
+     &       +(evdwij+eesij)*sssgrad*rrmij
           erij(1)=xj*rmij
           erij(2)=yj*rmij
           erij(3)=zj*rmij
@@ -4297,7 +4342,7 @@ cd        print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
 cd   &          (dcosg(k),k=1,3)
           do k=1,3
             ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*
-     &      fac_shield(i)**2*fac_shield(j)**2
+     &      fac_shield(i)**2*fac_shield(j)**2*sss
           enddo
 c          do k=1,3
 c            ghalf=0.5D0*ggg(k)
@@ -4317,11 +4362,11 @@ C                     print *,"before22", gelc_long(1,i), gelc_long(1,j)
           do k=1,3
             gelc(k,i)=gelc(k,i)
      &           +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
-     &           + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))
+     &           + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))*sss
      &           *fac_shield(i)**2*fac_shield(j)**2   
             gelc(k,j)=gelc(k,j)
      &           +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
-     &           + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))
+     &           + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))*sss
      &           *fac_shield(i)**2*fac_shield(j)**2
             gelc_long(k,j)=gelc_long(k,j)+ggg(k)
             gelc_long(k,i)=gelc_long(k,i)-ggg(k)
@@ -4557,7 +4602,7 @@ C           fac_shield(i)=0.4
 C           fac_shield(j)=0.6
           endif
           eel_loc_ij=eel_loc_ij
-     &    *fac_shield(i)*fac_shield(j)
+     &    *fac_shield(i)*fac_shield(j)*sss
 c          if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
 c     &            'eelloc',i,j,eel_loc_ij
 C Now derivative over eel_loc
@@ -4615,7 +4660,7 @@ C Calculate patrial derivative for theta angle
      &     +a23*gmuij1(2)
      &     +a32*gmuij1(3)
      &     +a33*gmuij1(4))
-     &    *fac_shield(i)*fac_shield(j)
+     &    *fac_shield(i)*fac_shield(j)*sss
 c         write(iout,*) "derivative over thatai"
 c         write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
 c     &   a33*gmuij1(4) 
@@ -4631,7 +4676,7 @@ c     &   a33*gmuij2(4)
      &     +a33*gmuij2(4)
          gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
      &      geel_loc_ij*wel_loc
-     &    *fac_shield(i)*fac_shield(j)
+     &    *fac_shield(i)*fac_shield(j)*sss
 
 c  Derivative over j residue
          geel_loc_ji=a22*gmuji1(1)
@@ -4644,7 +4689,7 @@ c     &   a33*gmuji1(4)
 
         gloc(nphi+j,icg)=gloc(nphi+j,icg)+
      &      geel_loc_ji*wel_loc
-     &    *fac_shield(i)*fac_shield(j)
+     &    *fac_shield(i)*fac_shield(j)*sss
 
          geel_loc_ji=
      &     +a22*gmuji2(1)
@@ -4656,7 +4701,7 @@ c         write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3),
 c     &   a33*gmuji2(4)
          gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+
      &      geel_loc_ji*wel_loc
-     &    *fac_shield(i)*fac_shield(j)
+     &    *fac_shield(i)*fac_shield(j)*sss
 #endif
 cd          write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
 
@@ -4672,17 +4717,21 @@ C Partial derivatives in virtual-bond dihedral angles gamma
      &    gel_loc_loc(i-1)=gel_loc_loc(i-1)+ 
      &            (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
      &           +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j))
-     &    *fac_shield(i)*fac_shield(j)
+     &    *fac_shield(i)*fac_shield(j)*sss
 
           gel_loc_loc(j-1)=gel_loc_loc(j-1)+ 
      &           (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
      &           +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j))
-     &    *fac_shield(i)*fac_shield(j)
+     &    *fac_shield(i)*fac_shield(j)*sss
 C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
+          aux=eel_loc_ij/sss*sssgrad*rmij
+          ggg(1)=aux*xj
+          ggg(2)=aux*yj
+          ggg(3)=aux*zj
           do l=1,3
-            ggg(l)=(agg(l,1)*muij(1)+
+            ggg(l)=ggg(l)+(agg(l,1)*muij(1)+
      &          agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4))
-     &    *fac_shield(i)*fac_shield(j)
+     &    *fac_shield(i)*fac_shield(j)*sss
             gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
             gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
 cgrad            ghalf=0.5d0*ggg(l)
@@ -4698,19 +4747,19 @@ C Remaining derivatives of eello
           do l=1,3
             gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+
      &        aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4))
-     &    *fac_shield(i)*fac_shield(j)
+     &    *fac_shield(i)*fac_shield(j)*sss
 
             gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+
      &     aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4))
-     &    *fac_shield(i)*fac_shield(j)
+     &    *fac_shield(i)*fac_shield(j)*sss
 
             gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+
      &       aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4))
-     &    *fac_shield(i)*fac_shield(j)
+     &    *fac_shield(i)*fac_shield(j)*sss
 
             gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+
      &     aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4))
-     &    *fac_shield(i)*fac_shield(j)
+     &    *fac_shield(i)*fac_shield(j)*sss
 
           enddo
           ENDIF
@@ -4800,9 +4849,9 @@ C                fac_shield(i)=0.4d0
 C                fac_shield(j)=0.6d0
                 endif
                 ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
-     &          *fac_shield(i)*fac_shield(j) 
+     &          *fac_shield(i)*fac_shield(j)*sss
                 ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
-     &          *fac_shield(i)*fac_shield(j)
+     &          *fac_shield(i)*fac_shield(j)*sss
 C Diagnostics. Comment out or remove after debugging!
 c               ees0p(num_conti,i)=0.5D0*fac3*ees0pij
 c               ees0m(num_conti,i)=0.5D0*fac3*ees0mij
@@ -4851,11 +4900,17 @@ cd              fprimcont=0.0D0
                   gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
                 enddo
                 gggp(1)=gggp(1)+ees0pijp*xj
+     &          +ees0p(num_conti,i)/sss*rmij*xj*sssgrad                
                 gggp(2)=gggp(2)+ees0pijp*yj
+     &          +ees0p(num_conti,i)/sss*rmij*yj*sssgrad
                 gggp(3)=gggp(3)+ees0pijp*zj
+     &          +ees0p(num_conti,i)/sss*rmij*zj*sssgrad
                 gggm(1)=gggm(1)+ees0mijp*xj
+     &          +ees0m(num_conti,i)/sss*rmij*xj*sssgrad                
                 gggm(2)=gggm(2)+ees0mijp*yj
+     &          +ees0m(num_conti,i)/sss*rmij*yj*sssgrad
                 gggm(3)=gggm(3)+ees0mijp*zj
+     &          +ees0m(num_conti,i)/sss*rmij*zj*sssgrad
 C Derivatives due to the contact function
                 gacont_hbr(1,num_conti,i)=fprimcont*xj
                 gacont_hbr(2,num_conti,i)=fprimcont*yj
@@ -4870,28 +4925,28 @@ cgrad                  ghalfm=0.5D0*gggm(k)
                   gacontp_hb1(k,num_conti,i)=!ghalfp
      &              +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
      &              + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
-     &          *fac_shield(i)*fac_shield(j)
+     &          *fac_shield(i)*fac_shield(j)*sss
 
                   gacontp_hb2(k,num_conti,i)=!ghalfp
      &              +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
      &              + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
-     &          *fac_shield(i)*fac_shield(j)
+     &          *fac_shield(i)*fac_shield(j)*sss
 
                   gacontp_hb3(k,num_conti,i)=gggp(k)
-     &          *fac_shield(i)*fac_shield(j)
+     &          *fac_shield(i)*fac_shield(j)*sss
 
                   gacontm_hb1(k,num_conti,i)=!ghalfm
      &              +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
      &              + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
-     &          *fac_shield(i)*fac_shield(j)
+     &          *fac_shield(i)*fac_shield(j)*sss
 
                   gacontm_hb2(k,num_conti,i)=!ghalfm
      &              +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
      &              + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
-     &          *fac_shield(i)*fac_shield(j)
+     &          *fac_shield(i)*fac_shield(j)*sss
 
                   gacontm_hb3(k,num_conti,i)=gggm(k)
-     &          *fac_shield(i)*fac_shield(j)
+     &          *fac_shield(i)*fac_shield(j)*sss
 
                 enddo
 C Diagnostics. Comment out or remove after debugging!
@@ -5672,7 +5727,7 @@ C This subroutine calculates the excluded-volume interaction energy between
 C peptide-group centers and side chains and its gradient in virtual-bond and
 C side-chain vectors.
 C
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
       include 'COMMON.GEO'
       include 'COMMON.VAR'
@@ -5685,7 +5740,14 @@ C
       include 'COMMON.CONTROL'
       include 'COMMON.SPLITELE'
       integer xshift,yshift,zshift
-      dimension ggg(3)
+      double precision ggg(3)
+      integer i,iint,j,k,iteli,itypj,subchap
+      double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1,
+     & fac,e1,e2,rij
+      double precision evdw2,evdw2_14,evdwij
+      double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,
+     & dist_temp, dist_init
+      double precision sscale,sscagrad
       evdw2=0.0D0
       evdw2_14=0.0d0
 c        print *,boxxsize,boxysize,boxzsize,'wymiary pudla'
@@ -5831,8 +5893,9 @@ c          if (sss.eq.0) print *,'czasem jest OK'
           endif
           evdwij=e1+e2
           evdw2=evdw2+evdwij*sss
-          if (energy_dec) write (iout,'(a6,2i5,0pf7.3,2i3,3e11.3)')
-     &        'evdw2',i,j,evdwij,iteli,itypj,fac,aad(itypj,iteli),
+          if (energy_dec) write (iout,'(a6,2i5,3f7.3,2i3,3e11.3)')
+     &        'evdw2',i,j,1.0d0/dsqrt(rrij),sss,
+     &       evdwij,iteli,itypj,fac,aad(itypj,iteli),
      &       bad(itypj,iteli)
 C
 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
index c73426c..6806e62 100644 (file)
@@ -46,6 +46,7 @@ cMS$ATTRIBUTES C ::  proc_proc
       include 'COMMON.DERIV'
       include 'COMMON.SPLITELE'
       include 'COMMON.VAR'
+      include 'COMMON.MD'
 c Common blocks from the diagonalization routines
       integer IR,IW,IP,IJK,IPK,IDAF,NAV,IODA,KDIAG,ICORFL,IXDR
       integer i,idumm,j,k,l,ichir1,ichir2,iblock,m
@@ -68,7 +69,7 @@ c NaNQ initialization
       call proc_proc(rr,i)
 #endif
 #endif
-
+      itime_mat=0.
       kdiag=0
       icorfl=0
       iw=2
index 4230e10..1180645 100644 (file)
@@ -891,6 +891,7 @@ c---------------------------------------------------------------------------
       include 'COMMON.IOUNITS'
       include 'COMMON.LAGRANGE.5diag'
       include 'COMMON.INTERACT'
+      include 'COMMON.VAR'
       integer ndim
       double precision forces(3*ndim),accel(3,0:maxres2),rs(ndim),
      &  xsolv(ndim),d_a_vec(6*nres)
@@ -904,9 +905,18 @@ Compute accelerations in Calpha and SC
           innt=chain_border(1,ichain)
           inct=chain_border(2,ichain)
           do i=iposc,iposc+n-1
-            rs(i)=forces(3*(i-1)+j)
+            rs(i-iposc+1)=forces(3*(i-1)+j)
           enddo
+#ifdef DEBUG
+          write (iout,*) "j",j," chain",ichain
+          write (iout,*) "rs"
+          write (iout,'(f10.5)') (rs(i),i=1,n)
+#endif
           call FDISYS (n,DM(iposc),DU1(iposc),DU2(iposc),rs,xsolv)
+#ifdef DEBUG
+          write (iout,*) "xsolv"
+          write (iout,'(f10.5)') (xsolv(i),i=1,n)
+#endif
           ind=1
           do i=innt,inct
             if (itype(i).eq.10)then
@@ -988,7 +998,7 @@ C Convert d_a to virtual-bon-vector basis
       enddo
 #ifdef DEBUG
       write (iout,*) "d_a_vec"
-      write (iout,'(3f10.5)') (d_a_vec(j),j=1,dimen3)
+      write (iout,'(3f10.5)') (d_a_vec(j),j=1,3*(nct-nnt+nside))
 #endif
       return
       end
diff --git a/source/unres/src-HCD-5D/make_xx_list.F b/source/unres/src-HCD-5D/make_xx_list.F
new file mode 100644 (file)
index 0000000..fb6c055
--- /dev/null
@@ -0,0 +1,489 @@
+      subroutine make_SCSC_inter_list
+      implicit none
+      include "DIMENSIONS"
+#ifdef MPI
+      include 'mpif.h'
+      include "COMMON.SETUP"
+#endif
+      include "COMMON.CHAIN"
+      include "COMMON.INTERACT"
+      include "COMMON.SPLITELE"
+      include "COMMON.IOUNITS"
+      double precision xi,yi,zi,xj,yj,zj,xj_safe,yj_safe,zj_safe,
+     &  xj_temp,yj_temp,zj_temp
+      double precision dist_init, dist_temp,r_buff_list
+      integer contlisti(200*maxres),contlistj(200*maxres)
+!      integer :: newcontlisti(200*nres),newcontlistj(200*nres) 
+      integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint,
+     &  ilist_sc,g_ilist_sc
+      integer displ(0:max_fg_procs),i_ilist_sc(0:max_fg_procs),ierr
+!            print *,"START make_SC"
+#ifdef DEBUG
+      write (iout,*) "make_SCSC_inter_list"
+#endif
+          r_buff_list=5.0d0
+            ilist_sc=0
+            do i=iatsc_s,iatsc_e
+             itypi=iabs(itype(i))
+             if (itypi.eq.ntyp1) cycle
+             xi=c(1,nres+i)
+             yi=c(2,nres+i)
+             zi=c(3,nres+i)
+             xi=dmod(xi,boxxsize)
+             if (xi.lt.0) xi=xi+boxxsize
+             yi=dmod(yi,boxysize)
+             if (yi.lt.0) yi=yi+boxysize
+             zi=dmod(zi,boxzsize)
+             if (zi.lt.0) zi=zi+boxzsize
+             do iint=1,nint_gr(i)
+              do j=istart(i,iint),iend(i,iint)
+               itypj=iabs(itype(j))
+               if (itypj.eq.ntyp1) cycle
+               xj=c(1,nres+j)
+               yj=c(2,nres+j)
+               zj=c(3,nres+j)
+               xj=dmod(xj,boxxsize)
+               if (xj.lt.0) xj=xj+boxxsize
+               yj=dmod(yj,boxysize)
+               if (yj.lt.0) yj=yj+boxysize
+               zj=dmod(zj,boxzsize)
+               if (zj.lt.0) zj=zj+boxzsize
+               dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+               xj_safe=xj
+               yj_safe=yj
+               zj_safe=zj
+               subchap=0
+               do xshift=-1,1
+               do yshift=-1,1
+               do zshift=-1,1
+               xj=xj_safe+xshift*boxxsize
+               yj=yj_safe+yshift*boxysize
+               zj=zj_safe+zshift*boxzsize
+               dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+               if(dist_temp.lt.dist_init) then
+                dist_init=dist_temp
+                xj_temp=xj
+                yj_temp=yj
+                zj_temp=zj
+                subchap=1
+               endif
+               enddo
+               enddo
+               enddo
+               if (subchap.eq.1) then
+               xj=xj_temp-xi
+               yj=yj_temp-yi
+               zj=zj_temp-zi
+               else
+               xj=xj_safe-xi
+               yj=yj_safe-yi
+               zj=zj_safe-zi
+               endif
+! r_buff_list is a read value for a buffer 
+               if (sqrt(dist_init).le.(r_cut_int+r_buff_list)) then
+! Here the list is created
+                 ilist_sc=ilist_sc+1
+! this can be substituted by cantor and anti-cantor
+                 contlisti(ilist_sc)=i
+                 contlistj(ilist_sc)=j
+
+               endif
+             enddo
+             enddo
+             enddo
+!         call MPI_Reduce(ilist_sc,g_ilist_sc,1,&
+!          MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
+!        call MPI_Gather(newnss,1,MPI_INTEGER,&
+!                        i_newnss,1,MPI_INTEGER,king,FG_COMM,IERR)
+#ifdef MPI
+#ifdef DEBUG
+      write (iout,*) "before MPIREDUCE",ilist_sc
+      do i=1,ilist_sc
+      write (iout,*) i,contlisti(i),contlistj(i)
+      enddo
+#endif
+      if (nfgtasks.gt.1)then
+
+        call MPI_Reduce(ilist_sc,g_ilist_sc,1,
+     &                  MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
+c        write(iout,*) "before bcast",g_ilist_sc
+        call MPI_Gather(ilist_sc,1,MPI_INTEGER,
+     &                  i_ilist_sc,1,MPI_INTEGER,king,FG_COMM,IERR)
+        displ(0)=0
+        do i=1,nfgtasks-1,1
+          displ(i)=i_ilist_sc(i-1)+displ(i-1)
+        enddo
+!        write(iout,*) "before gather",displ(0),displ(1)        
+        call MPI_Gatherv(contlisti,ilist_sc,MPI_INTEGER,
+     &                   newcontlisti,i_ilist_sc,displ,MPI_INTEGER,
+     &                   king,FG_COMM,IERR)
+        call MPI_Gatherv(contlistj,ilist_sc,MPI_INTEGER,
+     &                   newcontlistj,i_ilist_sc,displ,MPI_INTEGER,
+     &                   king,FG_COMM,IERR)
+        call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM,IERR)
+!        write(iout,*) "before bcast",g_ilist_sc
+!        call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+        call MPI_Bcast(newcontlisti,g_ilist_sc,MPI_INT,king,FG_COMM,
+     &                 IERR)
+        call MPI_Bcast(newcontlistj,g_ilist_sc,MPI_INT,king,FG_COMM,
+     &                 IERR)
+!        call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+
+        else
+#endif
+        g_ilist_sc=ilist_sc
+
+        do i=1,ilist_sc
+        newcontlisti(i)=contlisti(i)
+        newcontlistj(i)=contlistj(i)
+        enddo
+#ifdef MPI
+        endif
+#endif      
+#ifdef DEBUG
+      write (iout,*) "after GATHERV",g_ilist_sc
+      do i=1,g_ilist_sc
+      write (iout,*) i,newcontlisti(i),newcontlistj(i)
+      enddo
+#endif
+        call int_bounds(g_ilist_sc,g_listscsc_start,g_listscsc_end)
+      return
+      end
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+      subroutine make_SCp_inter_list
+      implicit none
+      include "DIMENSIONS"
+#ifdef MPI
+      include 'mpif.h'
+      include "COMMON.SETUP"
+#endif
+      include "COMMON.CHAIN"
+      include "COMMON.INTERACT"
+      include "COMMON.SPLITELE"
+      include "COMMON.IOUNITS"
+      double precision xi,yi,zi,xj,yj,zj,xj_safe,yj_safe,zj_safe,
+     &  xj_temp,yj_temp,zj_temp
+      double precision dist_init, dist_temp,r_buff_list
+      integer contlistscpi(200*maxres),contlistscpj(200*maxres)
+!      integer :: newcontlistscpi(200*nres),newcontlistscpj(200*nres)
+      integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint,
+     &  ilist_scp,g_ilist_scp
+      integer displ(0:max_fg_procs),i_ilist_scp(0:max_fg_procs),ierr
+      integer contlistscpi_f(200*maxres),contlistscpj_f(200*maxres)
+      integer ilist_scp_first,ifirstrun,g_ilist_sc
+!            print *,"START make_SC"
+#ifdef DEBUG
+      write (iout,*) "make_SCp_inter_list"
+#endif
+      r_buff_list=5.0
+            ilist_scp=0
+            ilist_scp_first=0
+      do i=iatscp_s,iatscp_e
+        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+        xi=0.5D0*(c(1,i)+c(1,i+1))
+        yi=0.5D0*(c(2,i)+c(2,i+1))
+        zi=0.5D0*(c(3,i)+c(3,i+1))
+          xi=mod(xi,boxxsize)
+          if (xi.lt.0) xi=xi+boxxsize
+          yi=mod(yi,boxysize)
+          if (yi.lt.0) yi=yi+boxysize
+          zi=mod(zi,boxzsize)
+          if (zi.lt.0) zi=zi+boxzsize
+
+        do iint=1,nscp_gr(i)
+
+        do j=iscpstart(i,iint),iscpend(i,iint)
+          itypj=iabs(itype(j))
+          if (itypj.eq.ntyp1) cycle
+! Uncomment following three lines for SC-p interactions
+!         xj=c(1,nres+j)-xi
+!         yj=c(2,nres+j)-yi
+!         zj=c(3,nres+j)-zi
+! Uncomment following three lines for Ca-p interactions
+!          xj=c(1,j)-xi
+!          yj=c(2,j)-yi
+!          zj=c(3,j)-zi
+          xj=c(1,j)
+          yj=c(2,j)
+          zj=c(3,j)
+          xj=mod(xj,boxxsize)
+          if (xj.lt.0) xj=xj+boxxsize
+          yj=mod(yj,boxysize)
+          if (yj.lt.0) yj=yj+boxysize
+          zj=mod(zj,boxzsize)
+          if (zj.lt.0) zj=zj+boxzsize
+      dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+      xj_safe=xj
+      yj_safe=yj
+      zj_safe=zj
+      subchap=0
+      do xshift=-1,1
+      do yshift=-1,1
+      do zshift=-1,1
+          xj=xj_safe+xshift*boxxsize
+          yj=yj_safe+yshift*boxysize
+          zj=zj_safe+zshift*boxzsize
+          dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+          if(dist_temp.lt.dist_init) then
+            dist_init=dist_temp
+            xj_temp=xj
+            yj_temp=yj
+            zj_temp=zj
+            subchap=1
+          endif
+       enddo
+       enddo
+       enddo
+       if (subchap.eq.1) then
+          xj=xj_temp-xi
+          yj=yj_temp-yi
+          zj=zj_temp-zi
+       else
+          xj=xj_safe-xi
+          yj=yj_safe-yi
+          zj=zj_safe-zi
+       endif
+#ifdef DEBUG
+                ! r_buff_list is a read value for a buffer 
+              if ((sqrt(dist_init).le.(r_cut_int)).and.(ifirstrun.eq.0))
+     &        then
+! Here the list is created
+                 ilist_scp_first=ilist_scp_first+1
+! this can be substituted by cantor and anti-cantor
+                 contlistscpi_f(ilist_scp_first)=i
+                 contlistscpj_f(ilist_scp_first)=j
+              endif
+#endif
+! r_buff_list is a read value for a buffer 
+               if (sqrt(dist_init).le.(r_cut_int+r_buff_list)) then
+! Here the list is created
+                 ilist_scp=ilist_scp+1
+! this can be substituted by cantor and anti-cantor
+                 contlistscpi(ilist_scp)=i
+                 contlistscpj(ilist_scp)=j
+              endif
+             enddo
+             enddo
+             enddo
+#ifdef MPI
+#ifdef DEBUG
+      write (iout,*) "before MPIREDUCE",ilist_scp
+      do i=1,ilist_scp
+      write (iout,*) i,contlistscpi(i),contlistscpj(i)
+      enddo
+#endif
+      if (nfgtasks.gt.1)then
+
+        call MPI_Reduce(ilist_scp,g_ilist_scp,1,
+     &    MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
+c        write(iout,*) "before bcast",g_ilist_sc
+        call MPI_Gather(ilist_scp,1,MPI_INTEGER,
+     &                  i_ilist_scp,1,MPI_INTEGER,king,FG_COMM,IERR)
+        displ(0)=0
+        do i=1,nfgtasks-1,1
+          displ(i)=i_ilist_scp(i-1)+displ(i-1)
+        enddo
+!        write(iout,*) "before gather",displ(0),displ(1)
+        call MPI_Gatherv(contlistscpi,ilist_scp,MPI_INTEGER,
+     &                   newcontlistscpi,i_ilist_scp,displ,MPI_INTEGER,
+     &                   king,FG_COMM,IERR)
+        call MPI_Gatherv(contlistscpj,ilist_scp,MPI_INTEGER,
+     &                   newcontlistscpj,i_ilist_scp,displ,MPI_INTEGER,
+     &                   king,FG_COMM,IERR)
+        call MPI_Bcast(g_ilist_scp,1,MPI_INT,king,FG_COMM,IERR)
+!        write(iout,*) "before bcast",g_ilist_sc
+!        call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+        call MPI_Bcast(newcontlistscpi,g_ilist_scp,MPI_INT,king,FG_COMM,
+     &                   IERR)
+        call MPI_Bcast(newcontlistscpj,g_ilist_scp,MPI_INT,king,FG_COMM,
+     &                   IERR)
+!        call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+        else
+#endif
+        g_ilist_scp=ilist_scp
+
+        do i=1,ilist_scp
+        newcontlistscpi(i)=contlistscpi(i)
+        newcontlistscpj(i)=contlistscpj(i)
+        enddo
+#ifdef MPI
+        endif
+#endif
+#ifdef DEBUG
+      write (iout,*) "after MPIREDUCE",g_ilist_scp
+      do i=1,g_ilist_scp
+      write (iout,*) i,newcontlistscpi(i),newcontlistscpj(i)
+      enddo
+
+!      if (ifirstrun.eq.0) ifirstrun=1
+!      do i=1,ilist_scp_first
+!       do j=1,g_ilist_scp
+!        if ((newcontlistscpi(j).eq.contlistscpi_f(i)).and.&
+!         (newcontlistscpj(j).eq.contlistscpj_f(i))) go to 126
+!        enddo
+!       print *,itime_mat,"ERROR matrix needs updating"
+!       print *,contlistscpi_f(i),contlistscpj_f(i)
+!  126  continue
+!      enddo
+#endif
+        call int_bounds(g_ilist_scp,g_listscp_start,g_listscp_end)
+
+      return
+      end 
+!-----------------------------------------------------------------------------
+      subroutine make_pp_inter_list
+      implicit none
+      include "DIMENSIONS"
+#ifdef MPI
+      include 'mpif.h'
+      include "COMMON.SETUP"
+#endif
+      include "COMMON.CHAIN"
+      include "COMMON.INTERACT"
+      include "COMMON.SPLITELE"
+      include "COMMON.IOUNITS"
+      double precision xi,yi,zi,xj,yj,zj,xj_safe,yj_safe,zj_safe,
+     &  xj_temp,yj_temp,zj_temp
+      double precision xmedj,ymedj,zmedj
+      double precision dist_init, dist_temp,r_buff_list,dxi,dyi,dzi,
+     &  xmedi,ymedi,zmedi
+      double precision dx_normi,dy_normi,dz_normi,dxj,dyj,dzj,
+     &  dx_normj,dy_normj,dz_normj
+      integer contlistppi(200*maxres),contlistppj(200*maxres)
+!      integer :: newcontlistppi(200*nres),newcontlistppj(200*nres)
+      integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint,
+     &  ilist_pp,g_ilist_pp
+      integer displ(0:max_fg_procs),i_ilist_pp(0:max_fg_procs),ierr
+!            print *,"START make_SC"
+#ifdef DEBUG
+      write (iout,*) "make_pp_inter_list"
+#endif
+      ilist_pp=0
+      r_buff_list=5.0
+      do i=iatel_s,iatel_e
+        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+        dxi=dc(1,i)
+        dyi=dc(2,i)
+        dzi=dc(3,i)
+        dx_normi=dc_norm(1,i)
+        dy_normi=dc_norm(2,i)
+        dz_normi=dc_norm(3,i)
+        xmedi=c(1,i)+0.5d0*dxi
+        ymedi=c(2,i)+0.5d0*dyi
+        zmedi=c(3,i)+0.5d0*dzi
+          xmedi=dmod(xmedi,boxxsize)
+          if (xmedi.lt.0) xmedi=xmedi+boxxsize
+          ymedi=dmod(ymedi,boxysize)
+          if (ymedi.lt.0) ymedi=ymedi+boxysize
+          zmedi=dmod(zmedi,boxzsize)
+          if (zmedi.lt.0) zmedi=zmedi+boxzsize
+             do j=ielstart(i),ielend(i)
+!          write (iout,*) i,j,itype(i),itype(j)
+          if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1) cycle
+! 1,j)
+          dxj=dc(1,j)
+          dyj=dc(2,j)
+          dzj=dc(3,j)
+          dx_normj=dc_norm(1,j)
+          dy_normj=dc_norm(2,j)
+          dz_normj=dc_norm(3,j)
+!          xj=c(1,j)+0.5D0*dxj-xmedi
+!          yj=c(2,j)+0.5D0*dyj-ymedi
+!          zj=c(3,j)+0.5D0*dzj-zmedi
+          xj=c(1,j)+0.5D0*dxj
+          yj=c(2,j)+0.5D0*dyj
+          zj=c(3,j)+0.5D0*dzj
+          xj=mod(xj,boxxsize)
+          if (xj.lt.0) xj=xj+boxxsize
+          yj=mod(yj,boxysize)
+          if (yj.lt.0) yj=yj+boxysize
+          zj=mod(zj,boxzsize)
+          if (zj.lt.0) zj=zj+boxzsize
+
+      dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
+      xj_safe=xj
+      yj_safe=yj
+      zj_safe=zj
+      do xshift=-1,1
+      do yshift=-1,1
+      do zshift=-1,1
+          xj=xj_safe+xshift*boxxsize
+          yj=yj_safe+yshift*boxysize
+          zj=zj_safe+zshift*boxzsize
+          dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
+          if(dist_temp.lt.dist_init) then
+            dist_init=dist_temp
+            xj_temp=xj
+            yj_temp=yj
+            zj_temp=zj
+          endif
+       enddo
+       enddo
+       enddo
+
+      if (sqrt(dist_init).le.(r_cut_int+r_buff_list)) then
+! Here the list is created
+                 ilist_pp=ilist_pp+1
+! this can be substituted by cantor and anti-cantor
+                 contlistppi(ilist_pp)=i
+                 contlistppj(ilist_pp)=j
+              endif
+             enddo
+             enddo
+!             enddo
+#ifdef MPI
+#ifdef DEBUG
+      write (iout,*) "before MPIREDUCE",ilist_pp
+      do i=1,ilist_pp
+      write (iout,*) i,contlistppi(i),contlistppj(i)
+      enddo
+#endif
+      if (nfgtasks.gt.1)then
+
+        call MPI_Reduce(ilist_pp,g_ilist_pp,1,
+     &    MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
+!        write(iout,*) "before bcast",g_ilist_sc
+        call MPI_Gather(ilist_pp,1,MPI_INTEGER,
+     &                  i_ilist_pp,1,MPI_INTEGER,king,FG_COMM,IERR)
+        displ(0)=0
+        do i=1,nfgtasks-1,1
+          displ(i)=i_ilist_pp(i-1)+displ(i-1)
+        enddo
+!        write(iout,*) "before gather",displ(0),displ(1)
+        call MPI_Gatherv(contlistppi,ilist_pp,MPI_INTEGER,
+     &                   newcontlistppi,i_ilist_pp,displ,MPI_INTEGER,
+     &                   king,FG_COMM,IERR)
+        call MPI_Gatherv(contlistppj,ilist_pp,MPI_INTEGER,
+     &                   newcontlistppj,i_ilist_pp,displ,MPI_INTEGER,
+     &                   king,FG_COMM,IERR)
+        call MPI_Bcast(g_ilist_pp,1,MPI_INT,king,FG_COMM,IERR)
+!        write(iout,*) "before bcast",g_ilist_sc
+!        call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+        call MPI_Bcast(newcontlistppi,g_ilist_pp,MPI_INT,king,FG_COMM,
+     &                   IERR)
+        call MPI_Bcast(newcontlistppj,g_ilist_pp,MPI_INT,king,FG_COMM,
+     &                   IERR)
+
+!        call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+
+        else
+#endif
+        g_ilist_pp=ilist_pp
+
+        do i=1,ilist_pp
+        newcontlistppi(i)=contlistppi(i)
+        newcontlistppj(i)=contlistppj(i)
+        enddo
+#ifdef MPI
+        endif
+#endif
+        call int_bounds(g_ilist_pp,g_listpp_start,g_listpp_end)
+#ifdef DEBUG
+      write (iout,*) "after MPIREDUCE",g_ilist_pp
+      do i=1,g_ilist_pp
+      write (iout,*) i,newcontlistppi(i),newcontlistppj(i)
+      enddo
+#endif
+      return
+      end
index 4765a41..8f4d05d 100644 (file)
@@ -171,6 +171,7 @@ c      call readi(controlcard,'IZ_SC',iz_sc,0)
       timlim=60.0D0*timlim
       safety = 60.0d0*safety
       modecalc=0
+      call readi(controlcard,"INTER_LIST_UPDATE",imatupdate,100)
       call reada(controlcard,"T_BATH",t_bath,300.0d0)
       minim=(index(controlcard,'MINIMIZE').gt.0)
       dccart=(index(controlcard,'CART').gt.0)
@@ -546,6 +547,10 @@ c  if performing umbrella sampling, fragments constrained are read from the frag
      &  "Initial time step of numerical integration:",d_time,
      &  " natural units"
        write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
+       write (iout,'(a60,f10.5,a)') "Cutoff on interactions",r_cut_int,
+     &   " A"
+       write(iout,'(a60,i5)')"Frequency of updating interaction list",
+     &   imatupdate
        if (RESPA) then
         write (iout,'(2a,i4,a)') 
      &    "A-MTS algorithm used; initial time step for fast-varying",
index b8069d9..09b6877 100644 (file)
@@ -70,6 +70,9 @@ c          write (iout,*) "friction_force j",j," ichain",ichain,
 c     &       " n",n," iposc",iposc,iposc+n-1
           innt=chain_border(1,ichain)
           inct=chain_border(2,ichain)
+c diagnostics
+c          innt=chain_border(1,1)
+c          inct=chain_border(2,1)
           do i=innt,inct
             vvec(ind+1)=v_work(j,i)
             ind=ind+1
@@ -324,6 +327,10 @@ c Compute the stochastic forces acting on bodies. Store in force.
         innt=chain_border(1,ichain)
         inct=chain_border(2,ichain)
         iposc=iposd_chain(ichain)
+c for debugging only
+c        innt=chain_border(1,1)
+c        inct=chain_border(2,1)
+c        iposc=iposd_chain(1)
 c        write (iout,*)"stochastic_force ichain=",ichain," innt",innt,
 c     &    " inct",inct," iposc",iposc
         do j=1,3