Adam's unres update
authorCezary Czaplewski <czarek@chem.univ.gda.pl>
Tue, 16 Feb 2021 14:11:53 +0000 (15:11 +0100)
committerCezary Czaplewski <czarek@chem.univ.gda.pl>
Tue, 16 Feb 2021 14:11:53 +0000 (15:11 +0100)
33 files changed:
source/unres/src-HCD-5D/COMMON.CHAIN
source/unres/src-HCD-5D/COMMON.CONTROL
source/unres/src-HCD-5D/COMMON.INTERACT
source/unres/src-HCD-5D/COMMON.MD
source/unres/src-HCD-5D/COMMON.TIME1
source/unres/src-HCD-5D/COMMON.TORCNSTR
source/unres/src-HCD-5D/DIMENSIONS
source/unres/src-HCD-5D/MD_A-MTS.F
source/unres/src-HCD-5D/cartprint.f
source/unres/src-HCD-5D/chain_symmetry.F
source/unres/src-HCD-5D/chainbuild.F
source/unres/src-HCD-5D/checkder_p.F
source/unres/src-HCD-5D/contact.f
source/unres/src-HCD-5D/dfa.F
source/unres/src-HCD-5D/elecont.f
source/unres/src-HCD-5D/energy_p_new-sep_barrier.F
source/unres/src-HCD-5D/energy_p_new_barrier.F
source/unres/src-HCD-5D/energy_split-sep.F
source/unres/src-HCD-5D/gen_rand_conf.F
source/unres/src-HCD-5D/gen_rand_conf_mchain.F [new file with mode: 0644]
source/unres/src-HCD-5D/geomout.F
source/unres/src-HCD-5D/gradient_p.F
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
source/unres/src-HCD-5D/minimize_p.F
source/unres/src-HCD-5D/orig_frame_chain.F [new file with mode: 0644]
source/unres/src-HCD-5D/readpdb-mult.F
source/unres/src-HCD-5D/readrtns_CSA.F
source/unres/src-HCD-5D/refsys.f
source/unres/src-HCD-5D/stochfric.F
source/unres/src-HCD-5D/timing.F
source/unres/src-HCD-5D/unres.F

index da83764..d666fd9 100644 (file)
@@ -1,7 +1,8 @@
       integer nres,nsup,nstart_sup,nz_start,nz_end,iz_sc,
      &  nres0,nstart_seq,nchain,chain_length,chain_border,iprzes,
-     &  chain_border1,ireschain,tabpermchain,npermchain,afmend,afmbeg,
-     &  nres_chomo,nmodel_start
+     &  chain_border1,ireschain,tabpermchain,npermchain,nequiv,
+     &  nchain_group,iequiv,mapchain,afmend,afmbeg,
+     &  nres_chomo,nmodel_start,nran_start
       double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm,t,r,
      & prod,rt,dc_work,cref,crefjlee,dc_norm2,velAFMconst,
      & totTafm,chomo
       common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2),
      & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2),
      & dc_norm2(3,0:maxres2),
-     & dc_work(MAXRES6),nres,nres0
+     & dc_work(MAXRES6),nres,nres0,nran_start
       common /rotmat/ t(3,3,maxres),r(3,3,maxres),prod(3,3,maxres),
      &                rt(3,3,maxres) 
       common /refstruct/ cref(3,maxres2+2),
      & crefjlee(3,maxres2+2),
      & nsup,nstart_sup,nstart_seq,iprzes,
      & chain_length(maxchain),npermchain,ireschain(maxres),
-     & tabpermchain(maxchain,maxperm),
+     & tabpermchain(maxchain,maxperm),nchain_group,
+     &  iequiv(maxchain,maxres),nequiv(maxchain),mapchain(maxchain),
      & chain_border(2,maxchain),chain_border1(2,maxchain),nchain
       common /from_zscore/ nz_start,nz_end,iz_sc
       double precision boxxsize,boxysize,boxzsize,enecut,sscut,
index da8581a..32ecc93 100644 (file)
@@ -2,7 +2,8 @@
 ! and output level.
 !... energy_dec = .true. means print energy decomposition matrix
       integer modecalc,iscode,indpdb,indback,indphi,iranconf,icheckgrad,
-     & inprint,i2ndstr,mucadyn,constr_dist,symetr,AFMlog,selfguide,
+     & inprint,i2ndstr,mucadyn,constr_dist,symetr,npermut,
+     & AFMlog,selfguide,
      & shield_mode,tor_mode,tubelog,constr_homology,homol_nset,
      & iprint
 !... minim = .true. means DO minimization.
@@ -20,6 +21,7 @@
      & unres_pdb,out_cart,out_int,vdisulf,searchsc,lmuca,dccart,mucadyn,
      & extconf,out1file,gmatout,selfguide,AFMlog,shield_mode,tor_mode,
      & tubelog,constr_dist,gnorm_check,gradout,split_ene,
-     & with_theta_constr,with_dihed_constr,symetr,usampl,loc_qlike,
+     & with_theta_constr,with_dihed_constr,symetr,npermut,usampl,loc_
+     & qlike,
      & adaptive,constr_homology,homol_nset,read2sigma,start_from_model,
      & read_homol_frag,out_template_coord,out_template_restr
index 14416ad..9b023e5 100644 (file)
      & 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
+C 10/26/20 Upgraded to interaction lists for RESPA
       integer newcontlisti(maxint_res*maxres),
+     & newcontlisti_long(maxint_res*maxres),
+     & newcontlisti_short(maxint_res*maxres),
      & newcontlistj(maxint_res*maxres),
+     & newcontlistj_long(maxint_res*maxres),
+     & newcontlistj_short(maxint_res*maxres),
      & newcontlistppi(maxint_res*maxres),
      & newcontlistppj(maxint_res*maxres),
-     & newcontlistpp_vdwi(maxint_res*maxres),
-     & newcontlistpp_vdwj(maxint_res*maxres),
+     & newcontlistpp_vdwi_short(maxint_res*maxres),
+     & newcontlistpp_vdwj_short(maxint_res*maxres),
      & newcontlistscpi(2*maxint_res*maxres),
+     & newcontlistscpi_long(2*maxint_res*maxres),
+     & newcontlistscpi_short(2*maxint_res*maxres),
      & newcontlistscpj(2*maxint_res*maxres),
+     & newcontlistscpj_long(2*maxint_res*maxres),
+     & newcontlistscpj_short(2*maxint_res*maxres),
      & g_listscsc_start,g_listscsc_end,g_listpp_start,g_listpp_end,
-     & g_listpp_vdw_start,g_listpp_vdw_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,newcontlistpp_vdwi,newcontlistpp_vdwj,
-     & g_listpp_vdw_start,g_listpp_vdw_end,
-     & newcontlistscpi,newcontlistscpj,g_listscp_start,
-     & g_listscp_end
+     & g_listscp_start,g_listscp_end,
+     & g_listscsc_start_long,g_listscsc_end_long,
+     & g_listscp_start_long,g_listscp_end_long,
+     & g_listscsc_start_short,g_listscsc_end_short,
+     & g_listpp_vdw_start_short,g_listpp_vdw_end_short,
+     & g_listscp_start_short,g_listscp_end_short
+      common /interact_list/
+     & newcontlisti,newcontlisti_long,newcontlisti_short,
+     & newcontlistj,newcontlistj_long,newcontlistj_short,
+     & g_listscsc_start,g_listscsc_start_long,g_listscsc_start_short,
+     & g_listscsc_end,g_listscsc_end_long,g_listscsc_end_short,
+     & newcontlistppi,newcontlistppj,g_listpp_start,
+     & g_listpp_end,
+     & newcontlistpp_vdwi_short,
+     & newcontlistpp_vdwj_short,
+     & g_listpp_vdw_start_short,
+     & g_listpp_vdw_end_short,
+     & newcontlistscpi,newcontlistscpi_long,newcontlistscpi_short,
+     & newcontlistscpj,newcontlistscpj_long,newcontlistscpj_short,
+     & g_listscp_start,g_listscp_start_long,g_listscp_start_short,
+     & g_listscp_end,g_listscp_end_long,g_listscp_end_short
 C 12/1/95 Array EPS included in the COMMON block.
       double precision eps,epslip,sigma,sigmaii,rs0,chi,chip,alp,
      & sigma0,sigii,
index cea18eb..c69e8e4 100644 (file)
@@ -3,13 +3,13 @@
      & dvmax,damax,edriftmax
       integer n_timestep,ntwx,ntwe,lang,count_reset_moment,
      & count_reset_vel,ntime_split,ntime_split0,
-     & maxtime_split,itime_mat,imatupdate
+     & maxtime_split,itime_mat,imatupdate,irest_freq
       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,itime_mat,imatupdate,
-     & ntwx,ntwe,lang,large,print_compon,tbf,rest,preminim,
+     & ntwx,ntwe,irest_freq,lang,large,print_compon,tbf,rest,preminim,
      & reset_moment,reset_vel,count_reset_moment,count_reset_vel,
      & rattle,RESPA
 ! Basic quantities 
index fc1e7d5..24d82ee 100644 (file)
@@ -15,7 +15,10 @@ c     FOUND_NAN - set by calcf to stop sumsl via stopx
      & time_vec,time_mat,time_ginvmult,time_fricmatmult,time_fric,
      & time_scatter_fmat,time_scatter_ginv,
      & time_fsample,time_scatter_fmatmult,time_scatter_ginvmult,
-     & time_stoch,t_eshort,t_elong,t_etotal,time_SAXS
+     & time_stoch,t_eshort,t_elong,t_etotal,time_SAXS,time_list,
+     & time_evdw,time_eelec,time_escp,time_evdw_short,time_eelec_short,
+     & time_escp_short,time_evdw_long,time_eelec_long,time_escp_long,
+     & time_escpcalc,time_escpsetup
       common /timing/ t_init,t_MDsetup,t_langsetup,
      & t_MD,t_enegrad,t_sdsetup,time_bcast,time_reduce,time_gather,
      & time_sendrecv,time_scatter,time_barrier_e,time_barrier_g,
@@ -25,4 +28,7 @@ c     FOUND_NAN - set by calcf to stop sumsl via stopx
      & time_vec,time_mat,time_ginvmult,time_fricmatmult,time_fric,
      & time_fsample,time_scatter_fmatmult,time_scatter_ginvmult,
      & time_scatter_fmat,time_scatter_ginv,
-     & time_stoch,t_eshort,t_elong,t_etotal,time_SAXS
+     & time_stoch,t_eshort,t_elong,t_etotal,time_SAXS,time_list,
+     & time_evdw,time_eelec,time_escp,time_evdw_short,time_eelec_short,
+     & time_escp_short,time_evdw_long,time_eelec_long,time_escp_long,
+     & time_escpcalc,time_escpsetup
index 9476b50..5219a92 100644 (file)
@@ -1,4 +1,5 @@
-      integer ndih_constr,idih_constr(maxdih_constr),ntheta_constr,
+      integer ndih_constr,idih_constr(maxdih_constr),
+     & iconstr_dih(maxres),ntheta_constr,
      & itheta_constr(maxdih_constr)
       integer ndih_nconstr,idih_nconstr(maxdih_constr)
       integer idihconstr_start,idihconstr_end,ithetaconstr_start,
@@ -10,7 +11,7 @@
      & vpsipred(3,maxdih_constr),sdihed(2,maxdih_constr),wdihc
       common /torcnstr/ phi0,drange,ftors,theta_constr0,theta_drange,
      & for_thet_constr,vpsipred,sdihed,wdihc,
-     &  ndih_constr,idih_constr,
+     &  ndih_constr,idih_constr,iconstr_dih,
      &  ndih_nconstr,idih_nconstr,idihconstr_start,idihconstr_end,
      & ntheta_constr,itheta_constr,ithetaconstr_start,
      & ithetaconstr_end,raw_psipred
index 4236776..3c2c924 100644 (file)
@@ -16,7 +16,8 @@ C Max. number of coarse-grain processors
       parameter (max_cg_procs=maxprocs)
 C Max. number of AA residues
       integer maxres
-      parameter (maxres=20000)
+c      parameter (maxres=70000)
+      parameter (maxres=1000)
 C Max. number of AA residues per chain
       integer maxres_chain
       parameter (maxres_chain=1200)
@@ -30,9 +31,10 @@ C Appr. max. number of interaction sites
      &           mmaxres2_chain=maxres2_chain*(maxres2_chain+1)/2)
 C Max number of symetric chains
       integer maxchain
-      parameter (maxchain=50)
+      parameter (maxchain=250)
       integer maxperm
-      parameter (maxperm=5040) 
+c      parameter (maxperm=5040) 
+      parameter (maxperm=3628800) 
 C Max. number of variables
       integer maxvar
       parameter (maxvar=6*maxres)
@@ -42,8 +44,8 @@ C Max. number of groups of interactions that a given SC is involved in
 C Max. number of derivatives of virtual-bond and side-chain vectors in theta
 C or phi.
       integer maxdim
-c      parameter (maxdim=(maxres_chain-1)*(maxres_chain-2)/2)
-      parameter (maxdim=(maxres-1)*(maxres-2)/2)
+      parameter (maxdim=(maxres_chain-1)*(maxres_chain-2)/2)
+c      parameter (maxdim=(maxres-1)*(maxres-2)/2)
 C Max. number of SC contacts
       integer maxcont
       parameter (maxcont=12*maxres)
@@ -60,8 +62,7 @@ C     include in template-based/contact distance restraints.
       parameter (maxcont_res=200)
 C Max. number of distance/contact-distance restraints
       integer maxdim_cont
-c      parameter (maxdim_cont=maxres*maxcont_res)
-      parameter (maxdim_cont=maxres*1000)
+      parameter (maxdim_cont=maxres*maxcont_res)
 C Number of AA types (at present only natural AA's will be handled
       integer ntyp,ntyp1
       parameter (ntyp=24,ntyp1=ntyp+1)
index e504cbd..8cebc8f 100644 (file)
@@ -65,6 +65,7 @@ c
       t_enegrad=0.0d0
       t_sdsetup=0.0d0
       write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
+      write (iout,*) "Restart frequency:",irest_freq
 #ifdef MPI
       tt0=MPI_Wtime()
 #else
@@ -261,7 +262,7 @@ C          call check_ecartint
              call cartout(totT)
           endif
         endif
-        if (rstcount.eq.1000.or.itime.eq.n_timestep) then
+        if (rstcount.eq.irest_freq.or.itime.eq.n_timestep) then
            open(irest2,file=rest2name,status='unknown')
            write(irest2,*) totT,EK,potE,totE,t_bath
            do i=0,2*nres-1
@@ -290,7 +291,7 @@ C          call check_ecartint
      & '  End of MD calculation  '
 #ifdef TIMING_ENE
       write (iout,*) "time for etotal",t_etotal," elong",t_elong,
-     &  " eshort",t_eshort
+     &  " eshort",t_eshort," list",time_list
       write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,
      & " time_fricmatmult",time_fricmatmult," time_fsample ",
      & time_fsample
@@ -429,7 +430,7 @@ c Calculate energy and forces
         call zerograd
         call etotal(potEcomp)
 ! AL 4/17/17: Reduce the steps if NaNs occurred.
-        if (potEcomp(0).gt.0.99e20 .or. isnan(potEcomp(0)).gt.0) then
+        if (potEcomp(0).gt.0.99e20 .or. isnan(potEcomp(0))) then
           d_time=d_time/2
           cycle
         endif
@@ -1672,6 +1673,7 @@ c  Set up the initial conditions of a MD simulation
       double precision etot
       logical fail
       integer i_start_models(0:nodes-1)
+      double precision potEcomp_long(0:Max_Ene)
       write (iout,*) "init_MD INDPDB",indpdb
       d_time0=d_time
 c      write(iout,*) "d_time", d_time
@@ -1810,7 +1812,7 @@ c      rest2name = prefix(:ilen(prefix))//'.rst'
       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,'(i7,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
      &   (d_t(j,i+nres),j=1,3)
        enddo
       endif
@@ -2080,6 +2082,7 @@ C 8/22/17 AL Minimize initial structure
       call zerograd
       call etotal(potEcomp)
       if (large) call enerprint(potEcomp)
+      if (RESPA) call etotal_long(potEcomp_long)
 #ifdef TIMING_ENE
 #ifdef MPI
       t_etotal=t_etotal+MPI_Wtime()-tt0
@@ -2092,7 +2095,7 @@ c      write (iout,*) "PotE-homology",potE-potEcomp(27)
       call cartgrad
       call lagrangian
       call max_accel
-      if (amax*d_time .gt. dvmax) then
+      if (amax*d_time .gt. dvmax .and. .not. respa) then
         d_time=d_time*dvmax/amax
         if(me.eq.king.or..not.out1file) write (iout,*) 
      &   "Time step reduced to",d_time,
index 339a89d..8a747b4 100644 (file)
@@ -15,6 +15,6 @@
      &          '     centroid coordinates'/
      1          '       ', 7X,'X',11X,'Y',11X,'Z',
      &                          10X,'X',11X,'Y',11X,'Z')
-  110 format (a,'(',i4,')',6f12.5)
+  110 format (a,'(',i5,')',6f12.5)
       return
       end  
index 8c36855..d10c33f 100644 (file)
@@ -1,5 +1,6 @@
       subroutine chain_symmetry(nchain,nres,itype,chain_border,
-     &    chain_length,npermchain,tabpermchain)
+     &    chain_length,npermchain,tabpermchain,nchain_group,nequiv,
+     &    iequiv,chaingroup)
 c
 c Determine chain symmetry. nperm is the number of permutations and
 c tabperchain contains the allowed permutations of the chains.
@@ -12,13 +13,17 @@ c
      &  chain_length(nchain),itemp(maxchain),
      &  npermchain,tabpermchain(maxchain,maxperm),
      &  tabperm(maxchain,maxperm),mapchain(maxchain),
-     &  iequiv(maxchain,maxres),iflag(maxres)
+     &  chaingroup(maxchain),iequiv(maxchain,maxres),iflag(maxres)
       integer i,j,k,l,ii,nchain_group,nequiv(maxchain),iieq,
      &  nperm,npermc,ind
+      logical lprn /.false./
       if (nchain.eq.1) then
         npermchain=1
         tabpermchain(1,1)=1
 c        print*,"npermchain",npermchain," tabpermchain",tabpermchain(1,1)
+        nchain_group=1
+        iequiv(1,1)=1
+        chaingroup(1)=1 
         return
       endif
 c
@@ -73,12 +78,24 @@ c            k=k+1
         do j=1,nequiv(i)
           ind=ind+1
           mapchain(ind)=iequiv(j,i)
+          chaingroup(ind)=i
         enddo
       enddo
       write (iout,*) "mapchain"
       do i=1,nchain
         write (iout,*) i,mapchain(i)
       enddo 
+      write (iout,*) "chaingroup"
+      do i=1,nchain
+        write (iout,*) i,chaingroup(i)
+      enddo 
+      if (npermut.eq.0) then
+        npermchain=1
+        do i=1,nchain
+          tabpermchain(i,1)=i
+        enddo
+        return
+      endif
       ii=0
       do i=1,nchain_group
         call permut(nequiv(i),nperm,tabperm)
@@ -117,10 +134,12 @@ c            k=k+1
         enddo
       enddo 
       write(iout,*) "Number of chain permutations",npermchain
+      if (lprn) then
       write(iout,*) "Permutations"
       do i=1,npermchain
         write(iout,'(20i4)') (tabpermchain(j,i),j=1,nchain)
       enddo
+      endif
       return
       end
 c---------------------------------------------------------------------
index 7902f15..a60e2bd 100644 (file)
@@ -108,11 +108,11 @@ C
       include 'COMMON.VAR'
       cost=dcos(theta(3))
       sint=dsin(theta(3))
-      t(1,1,1)=-cost
+      t(1,1,1)=cost
       t(1,2,1)=-sint 
       t(1,3,1)= 0.0D0
-      t(2,1,1)=-sint
-      t(2,2,1)= cost
+      t(2,1,1)=sint
+      t(2,2,1)=cost
       t(2,3,1)= 0.0D0
       t(3,1,1)= 0.0D0
       t(3,2,1)= 0.0D0
@@ -170,7 +170,7 @@ C
 C Locate CA(i) and SC(i-1)
 C
       implicit none
-      integer i,j
+      integer i,j,k
       double precision theti,phii,cost,sint,cosphi,sinphi
       include 'DIMENSIONS'
       include 'COMMON.CHAIN'
@@ -196,14 +196,15 @@ C
       sint=dsin(theti)
       cosphi=dcos(phii)
       sinphi=dsin(phii)
+c      write (iout,*) "locate_next_res i",i
 * Define the matrices of the rotation about the virtual-bond valence angles
 * theta, T(i,j,k), virtual-bond dihedral angles gamma (miscalled PHI in this
 * program), R(i,j,k), and, the cumulative matrices of rotation RT
       t(1,1,i-2)=-cost
-      t(1,2,i-2)=-sint 
+      t(1,2,i-2)= sint 
       t(1,3,i-2)= 0.0D0
       t(2,1,i-2)=-sint
-      t(2,2,i-2)= cost
+      t(2,2,i-2)=-cost
       t(2,3,i-2)= 0.0D0
       t(3,1,i-2)= 0.0D0
       t(3,2,i-2)= 0.0D0
@@ -212,26 +213,36 @@ C
       r(1,2,i-2)= 0.0D0
       r(1,3,i-2)= 0.0D0
       r(2,1,i-2)= 0.0D0
-      r(2,2,i-2)=-cosphi
-      r(2,3,i-2)= sinphi
+      r(2,2,i-2)= cosphi
+      r(2,3,i-2)=-sinphi
       r(3,1,i-2)= 0.0D0
       r(3,2,i-2)= sinphi
       r(3,3,i-2)= cosphi
       rt(1,1,i-2)=-cost
-      rt(1,2,i-2)=-sint
+      rt(1,2,i-2)= sint
       rt(1,3,i-2)=0.0D0
-      rt(2,1,i-2)=sint*cosphi
+      rt(2,1,i-2)=-sint*cosphi
       rt(2,2,i-2)=-cost*cosphi
-      rt(2,3,i-2)=sinphi
+      rt(2,3,i-2)=-sinphi
       rt(3,1,i-2)=-sint*sinphi
-      rt(3,2,i-2)=cost*sinphi
-      rt(3,3,i-2)=cosphi
+      rt(3,2,i-2)=-cost*sinphi
+      rt(3,3,i-2)= cosphi
       call matmult(prod(1,1,i-2),rt(1,1,i-2),prod(1,1,i-1))
+c      write (iout,*) "prod",i-2
+c      do j=1,3
+c        write (iout,*) (prod(j,k,i-2),k=1,3)
+c      enddo
+c      write (iout,*) "prod",i-1
+c      do j=1,3
+c        write (iout,*) (prod(j,k,i-1),k=1,3)
+c      enddo
       do j=1,3
         dc_norm(j,i-1)=prod(j,1,i-1)
         dc(j,i-1)=vbld(i)*prod(j,1,i-1)
         c(j,i)=c(j,i-1)+dc(j,i-1)
       enddo
+c      write (iout,*) "dc",i-1,(dc(j,i-1),j=1,3)
+c      write (iout,*) "c",i,(dc(j,i),j=1,3)
 cd    print '(2i3,2(3f10.5,5x))', i-1,i,(dc(j,i-1),j=1,3),(c(j,i),j=1,3)
 C 
 C Now calculate the coordinates of SC(i-1)
index e1448db..d0032da 100644 (file)
@@ -364,6 +364,7 @@ c-------------------------------------------------------------------------
       double precision dnorm1,dnorm2,be
       double precision time00
       double precision dist,alpha,beta
+      double precision time01
       if (lprn) write (iout,'(/a)') 'Recalculated internal coordinates'
 #ifdef TIMING
       time01=MPI_Wtime()
index f986380..ead0332 100644 (file)
@@ -137,7 +137,7 @@ c     print *,'nnt=',nnt,' nct=',nct
           i2=icont(2,i)
           it1=itype(i1)
           it2=itype(i2)
-          write (iout,'(i5,2x,a,i5,2x,a,i4)') 
+          write (iout,'(i10,2x,a,i7,2x,a,i7)') 
      &     i,restyp(it1),i1,restyp(it2),i2 
         enddo
       endif
@@ -184,8 +184,8 @@ c      enddo
         ii1=iharp(3,i)
         jj1=iharp(4,i)
         write (iout,*)
-        write (iout,'(20(a,i5,1x))') (restyp(itype(k)),k,k=i1,ii1)
-        write (iout,'(20(a,i5,1x))') (restyp(itype(k)),k,k=j1,jj1,-1)
+        write (iout,'(20(a,i7,1x))') (restyp(itype(k)),k,k=i1,ii1)
+        write (iout,'(20(a,i7,1x))') (restyp(itype(k)),k,k=j1,jj1,-1)
 c        do k=jj1,j1,-1
 c         write (iout,'(a,i3,$)') restyp(itype(k)),k
 c        enddo
index 1af0b44..3982b6d 100644 (file)
@@ -368,7 +368,6 @@ C     END OF BETA RESTRAINT
       edfadis=0
       gdfad=0.0d0
 
-c      write (2,*) "edfad",idfadis_start,idfadis_end
       do i=idfadis_start,idfadis_end
 
          iatm1=idislis(1,i)+ishiftca
index bf9056a..7c024ea 100644 (file)
@@ -9,7 +9,7 @@
       include 'COMMON.NAMES'
       logical lprint
       double precision elpp_6(2,2),elpp_3(2,2),ael6_(2,2),ael3_(2,2)
-      double precision app_(2,2),bpp_(2,2),rpp_(2,2)
+      double precision app_(2,2),bpp_(2,2),epp_(2,2),rpp_(2,2)
       integer ncont,icont(2,maxcont)
       double precision econt(maxcont)
       double precision boxshift
@@ -20,7 +20,7 @@
 *
 * as of 7/06/91.
 *
-c      data epp    / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/
+c      data epp_    / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/
       data rpp_    / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
       data elpp_6  /-0.2379d0,-0.2056d0,-0.2056d0,-0.0610d0/
       data elpp_3  / 0.0503d0, 0.0000d0, 0.0000d0, 0.0692d0/
@@ -55,7 +55,6 @@ c      data epp    / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/
         zmedi=zi+0.5*dzi
         call to_box(xmedi,ymedi,zmedi)
 c        write (iout,*) "i",xmedi,ymedi,zmedi
-c        write (iout,*) "i",xmedi,ymedi,zmedi
         do 4 j=i+2,nct-1
 c          write (iout,*) "i",i," j",j
           if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) goto 4
@@ -74,11 +73,12 @@ c          write (iout,*) "i",i," j",j
           xj=c(1,j)+0.5*dxj
           yj=c(2,j)+0.5*dyj
           zj=c(3,j)+0.5*dzj
-c          write (iout,*) "j",xj,yj,zj
           call to_box(xj,yj,zj)
-          xj=boxshift(xj-xi,boxxsize)
-          yj=boxshift(yj-yi,boxysize)
-          zj=boxshift(zj-zi,boxzsize)
+c          write (iout,*) "j",xj,yj,zj
+          xj=boxshift(xj-xmedi,boxxsize)
+          yj=boxshift(yj-ymedi,boxysize)
+          zj=boxshift(zj-zmedi,boxzsize)
+c          write (iout,*) "j",xj,yj,zj
           rij=xj*xj+yj*yj+zj*zj
           rrmij=1.0/(xj*xj+yj*yj+zj*zj)
           rmij=sqrt(rrmij)
@@ -105,7 +105,7 @@ c          write (iout,*) "j",xj,yj,zj
             econt(ncont)=eesij
           endif
           ees=ees+eesij
-c          write (iout,*) "i"," j",j," rij",dsqrt(rij)," eesij",eesij
+c          write (iout,*) "i",i," j",j," rij",dsqrt(rij)," eesij",eesij
     4   continue
     1 continue
       if (lprint) then
@@ -118,7 +118,7 @@ c          write (iout,*) "i"," j",j," rij",dsqrt(rij)," eesij",eesij
           i2=icont(2,i)
           it1=itype(i1)
           it2=itype(i2)
-          write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)')
+          write (iout,'(i7,2x,a,i6,2x,a,i6,f10.5)')
      &     i,restyp(it1),i1,restyp(it2),i2,econt(i)
         enddo
       endif
@@ -204,7 +204,7 @@ c              enddo
           i2=icont(2,i)
           it1=itype(i1)
           it2=itype(i2)
-          write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)')
+          write (iout,'(i7,2x,a,i6,2x,a,i6,f10.5)')
      &     i,restyp(it1),i1,restyp(it2),i2,econt(i)
         enddo
       endif
@@ -267,7 +267,7 @@ cd            write (iout,*) i1,j1,not_done
             ii1=max0(ii1-1,1)
             jj1=max0(jj1-1,1)
             nbeta=nbeta+1
-            if(lprint)write(iout,'(a,i3,4i4)')'parallel beta',
+            if(lprint)write(iout,'(a,i7,4i6)')'parallel beta',
      &               nbeta,ii1,i1,jj1,j1
 
             nbfrag=nbfrag+1
@@ -442,29 +442,29 @@ cd            write (iout,*) i1,j1,not_done
 
 
            if (lprint) then
-            write (iout,'(a,i3,4i4)')'antiparallel beta',
+            write (iout,'(a,i3,4i6)')'antiparallel beta',
      &                   nbeta,ii1-1,i1,jj1,j1-1
             nstrand=nstrand+1
             if (nstrand.le.9) then
-              write(12,'(a18,i1,a9,i3,a2,i3,a1)') 
+              write(12,'(a18,i1,a9,i6,a2,i6,a1)') 
      &          "DefPropRes 'strand",nstrand,
      &          "' 'num = ",ii1-2,"..",i1-1,"'"
             else
-              write(12,'(a18,i2,a9,i3,a2,i3,a1)') 
+              write(12,'(a18,i2,a9,i6,a2,i6,a1)') 
      &          "DefPropRes 'strand",nstrand,
      &          "' 'num = ",ii1-2,"..",i1-1,"'"
             endif
             nstrand=nstrand+1
             if (nstrand.le.9) then
-              write(12,'(a18,i1,a9,i3,a2,i3,a1)') 
+              write(12,'(a18,i1,a9,i6,a2,i6,a1)') 
      &          "DefPropRes 'strand",nstrand,
      &          "' 'num = ",j1-2,"..",jj1-1,"'"
             else
-              write(12,'(a18,i2,a9,i3,a2,i3,a1)') 
+              write(12,'(a18,i2,a9,i6,a2,i6,a1)') 
      &          "DefPropRes 'strand",nstrand,
      &          "' 'num = ",j1-2,"..",jj1-1,"'"
             endif
-              write(12,'(a8,4i4)')
+              write(12,'(a8,4i6)')
      &          "SetNeigh",ii1-2,i1-1,jj1-1,j1-2
            endif
           endif
index c4e54bc..f92aebb 100644 (file)
@@ -94,9 +94,9 @@ c      write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
       gg_lipi=0.0d0
       gg_lipj=0.0d0
 c      do i=iatsc_s,iatsc_e
-      do ikont=g_listscsc_start,g_listscsc_end
-        i=newcontlisti(ikont)
-        j=newcontlistj(ikont)
+      do ikont=g_listscsc_start_long,g_listscsc_end_long
+        i=newcontlisti_long(ikont)
+        j=newcontlistj_long(ikont)
         itypi=iabs(itype(i))
         if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
@@ -219,9 +219,9 @@ c      write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
       gg_lipi=0.0d0
       gg_lipj=0.0d0
 c      do i=iatsc_s,iatsc_e
-      do ikont=g_listscsc_start,g_listscsc_end
-        i=newcontlisti(ikont)
-        j=newcontlistj(ikont)
+      do ikont=g_listscsc_start_short,g_listscsc_end_short
+        i=newcontlisti_short(ikont)
+        j=newcontlistj_short(ikont)
         itypi=iabs(itype(i))
         if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
@@ -340,9 +340,9 @@ c     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
       gg_lipi=0.0d0
       gg_lipj=0.0d0
 c      do i=iatsc_s,iatsc_e
-      do ikont=g_listscsc_start,g_listscsc_end
-        i=newcontlisti(ikont)
-        j=newcontlistj(ikont)
+      do ikont=g_listscsc_start_long,g_listscsc_end_long
+        i=newcontlisti_long(ikont)
+        j=newcontlistj_long(ikont)
         itypi=iabs(itype(i))
         if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
@@ -462,9 +462,9 @@ c     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
       gg_lipi=0.0d0
       gg_lipj=0.0d0
 c      do i=iatsc_s,iatsc_e
-      do ikont=g_listscsc_start,g_listscsc_end
-        i=newcontlisti(ikont)
-        j=newcontlistj(ikont)
+      do ikont=g_listscsc_start_short,g_listscsc_end_short
+        i=newcontlisti_short(ikont)
+        j=newcontlistj_short(ikont)
         itypi=iabs(itype(i))
         if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
@@ -586,9 +586,9 @@ c     else
 c     endif
       ind=0
 c      do i=iatsc_s,iatsc_e
-      do ikont=g_listscsc_start,g_listscsc_end
-        i=newcontlisti(ikont)
-        j=newcontlistj(ikont)
+      do ikont=g_listscsc_start_long,g_listscsc_end_long
+        i=newcontlisti_long(ikont)
+        j=newcontlistj_long(ikont)
         itypi=iabs(itype(i))
         if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
@@ -734,9 +734,9 @@ c     else
 c     endif
       ind=0
 c      do i=iatsc_s,iatsc_e
-      do ikont=g_listscsc_start,g_listscsc_end
-        i=newcontlisti(ikont)
-        j=newcontlistj(ikont)
+      do ikont=g_listscsc_start_short,g_listscsc_end_short
+        i=newcontlisti_short(ikont)
+        j=newcontlistj_short(ikont)
         itypi=iabs(itype(i))
         if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
@@ -872,7 +872,7 @@ C
      & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
       double precision dist,sscale,sscagrad,sscagradlip,sscalelip
       double precision subchap,sss1,sssgrad1
-      double precision boxshift
+      double precision boxshift,rij1
       evdw=0.0D0
 ccccc      energy_dec=.false.
 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
@@ -882,9 +882,12 @@ c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
 c     if (icall.eq.0) lprn=.false.
       ind=0
 c      do i=iatsc_s,iatsc_e
-      do ikont=g_listscsc_start,g_listscsc_end
-        i=newcontlisti(ikont)
-        j=newcontlistj(ikont)
+      if (energy_dec)
+     & write(2,*) "g_listscsc_start_long,g_listscsc_end_long",
+     & g_listscsc_start_long,g_listscsc_end_long
+      do ikont=g_listscsc_start_long,g_listscsc_end_long
+        i=newcontlisti_long(ikont)
+        j=newcontlistj_long(ikont)
         itypi=iabs(itype(i))
         if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
@@ -940,9 +943,13 @@ c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
             dzj=dc_norm(3,nres+j)
             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
             rij=dsqrt(rrij)
-            sss1=sscale(1.0d0/rij,r_cut_int)
+            rij1=1.0d0/rij
+c            sss1=sscale(1.0d0/rij,r_cut_int)
+            sss1=sscale(rij1,r_cut_int)
             if (sss1.eq.0.0d0) cycle
-            sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
+            rij1=rij1/sigmaii(itypi,itypj)
+            sss=sscale(rij1,r_cut_respa)
+c            sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
             if (sss.lt.1.0d0) then
 C Calculate angle-dependent terms of energy and contributions to their
 C derivatives.
@@ -1026,6 +1033,7 @@ C This subroutine calculates the interaction energy of nonbonded side chains
 C assuming the Gay-Berne potential of interaction.
 C
       implicit none
+      include 'mpif.h'
       include 'DIMENSIONS'
       include 'COMMON.GEO'
       include 'COMMON.VAR'
@@ -1038,6 +1046,7 @@ C
       include 'COMMON.CALC'
       include 'COMMON.CONTROL'
       include "COMMON.SPLITELE"
+      include 'COMMON.TIME1'
       logical lprn
       double precision evdw
       integer itypi,itypj,itypi1,iint,ind,ikont
@@ -1046,6 +1055,8 @@ C
      & sslipj,ssgradlipj,ssgradlipi,sig,rij_shift,faclip
       double precision dist,sscale,sscagrad,sscagradlip,sscalelip
       double precision boxshift
+      double precision time01
+c      time01=MPI_Wtime()
       evdw=0.0D0
 ccccc      energy_dec=.false.
 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
@@ -1055,9 +1066,12 @@ c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
 c     if (icall.eq.0) lprn=.false.
       ind=0
 c      do i=iatsc_s,iatsc_e
-      do ikont=g_listscsc_start,g_listscsc_end
-        i=newcontlisti(ikont)
-        j=newcontlistj(ikont)
+      if (energy_dec)
+     & write(2,*) "g_listscsc_start_short,g_listscsc_end_short",
+     & g_listscsc_start_short,g_listscsc_end_short
+      do ikont=g_listscsc_start_short,g_listscsc_end_short
+        i=newcontlisti_short(ikont)
+        j=newcontlistj_short(ikont)
         itypi=iabs(itype(i))
         if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
@@ -1119,8 +1133,8 @@ c     &        (2.0d0-sslipi-sslipj)/2.0d0
             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
             rij=dsqrt(rrij)
             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
-          sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa)
             if (sss.gt.0.0d0) then
+          sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa)
 
 C Calculate angle-dependent terms of energy and contributions to their
 C derivatives.
@@ -1189,6 +1203,7 @@ C Calculate angular part of the gradient.
 c          enddo      ! j
 c        enddo        ! iint
       enddo          ! i
+c      time_evdw_short=time_evdw_short+MPI_Wtime()-time01
 c      write (iout,*) "Number of loop steps in EGB:",ind
 cccc      energy_dec=.false.
       return
@@ -1230,9 +1245,9 @@ c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
 c     if (icall.eq.0) lprn=.true.
       ind=0
 c      do i=iatsc_s,iatsc_e
-      do ikont=g_listscsc_start,g_listscsc_end
-        i=newcontlisti(ikont)
-        j=newcontlistj(ikont)
+      do ikont=g_listscsc_start_long,g_listscsc_end_long
+        i=newcontlisti_long(ikont)
+        j=newcontlistj_long(ikont)
         itypi=iabs(itype(i))
         if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
@@ -1392,9 +1407,9 @@ c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
 c     if (icall.eq.0) lprn=.true.
       ind=0
 c      do i=iatsc_s,iatsc_e
-      do ikont=g_listscsc_start,g_listscsc_end
-        i=newcontlisti(ikont)
-        j=newcontlistj(ikont)
+      do ikont=g_listscsc_start_short,g_listscsc_end_short
+        i=newcontlisti_short(ikont)
+        j=newcontlistj_short(ikont)
         itypi=iabs(itype(i))
         if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
@@ -1742,6 +1757,9 @@ c
 c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
 c
 c      do i=iatel_s,iatel_e
+      if (energy_dec)
+     & write(iout,*) "g_listpp_start,g_listpp_end",
+     & g_listpp_start,g_listpp_end
       do ikont=g_listpp_start,g_listpp_end
         i=newcontlistppi(ikont)
         j=newcontlistppj(ikont)
@@ -2812,9 +2830,12 @@ c      write (iout,*) "iatel_s_vdw",iatel_s_vdw,
 c     & " iatel_e_vdw",iatel_e_vdw
 c      call flush(iout)
 c      do i=iatel_s_vdw,iatel_e_vdw
-      do ikont=g_listpp_vdw_start,g_listpp_vdw_end
-        i=newcontlistpp_vdwi(ikont)
-        j=newcontlistpp_vdwj(ikont)
+      if (energy_dec)
+     & write(iout,*) "g_listpp_vdw_start_short,g_listpp_vdw_end_short",
+     & g_listpp_vdw_start_short,g_listpp_vdw_end_short
+      do ikont=g_listpp_vdw_start_short,g_listpp_vdw_end_short
+        i=newcontlistpp_vdwi_short(ikont)
+        j=newcontlistpp_vdwj_short(ikont)
         if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1) cycle
         dxi=dc(1,i)
         dyi=dc(2,i)
@@ -2939,9 +2960,12 @@ c      if (lprint_short)
 c     &  write (iout,*) 'ESCP_LONG iatscp_s=',iatscp_s,
 c     & ' iatscp_e=',iatscp_e
 c      do i=iatscp_s,iatscp_e
-      do ikont=g_listscp_start,g_listscp_end
-        i=newcontlistscpi(ikont)
-        j=newcontlistscpj(ikont)
+      if (energy_dec)
+     & write(iout,*)"g_listscp_start_long,g_listscp_end_long",
+     & g_listscp_start_long,g_listscp_end_long
+      do ikont=g_listscp_start_long,g_listscp_end_long
+        i=newcontlistscpi_long(ikont)
+        j=newcontlistscpj_long(ikont)
         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))
       double precision ggg(3)
       double precision sscale,sscagrad
       double precision boxshift
+      integer ikont
       evdw2=0.0D0
       evdw2_14=0.0d0
 cd    print '(a)','Enter ESCP'
 c      if (lprint_short) 
 c     &  write (iout,*) 'ESCP_SHORT iatscp_s=',iatscp_s,
 c     & ' iatscp_e=',iatscp_e
-      if (energy_dec) write (iout,*) "escp_short:",r_cut_int,rlamb
-      do i=iatscp_s,iatscp_e
+c      if (energy_dec) write (iout,*) "escp_short:",r_cut_int,rlamb
+      if (energy_dec)
+     & write(iout,*) "g_listscp_start_short,g_listscp_end_short",
+     & g_listscp_start_short,g_listscp_end_short
+      do ikont=g_listscp_start_short,g_listscp_end_short
+        i=newcontlistscpi_short(ikont)
+        j=newcontlistscpj_short(ikont)
         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))
@@ -3083,9 +3113,9 @@ c     & ' iatscp_e=',iatscp_e
 c        if (lprint_short) 
 c     &    write (iout,*) "i",i," itype",itype(i),itype(i+1),
 c     &     " nscp_gr",nscp_gr(i)   
-        do iint=1,nscp_gr(i)
-
-        do j=iscpstart(i,iint),iscpend(i,iint)
+c        do iint=1,nscp_gr(i)
+c
+c        do j=iscpstart(i,iint),iscpend(i,iint)
           itypj=iabs(itype(j))
 c        if (lprint_short)
 c     &    write (iout,*) "j",j," itypj",itypj
@@ -3149,9 +3179,9 @@ c             gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
               gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
             enddo
           endif
-        enddo
+c        enddo
 
-        enddo ! iint
+c        enddo ! iint
       enddo ! i
       do i=1,nct
         do j=1,3
index 6d6a817..2d94dc0 100644 (file)
@@ -37,6 +37,7 @@ c      include 'COMMON.MD'
      & eliptran,Eafmforce,Etube,
      & esaxs_constr,ehomology_constr,edfator,edfanei,edfabet
       integer n_corr,n_corr1
+      double precision time01
 #ifdef MPI      
 c      print*,"ETOTAL Processor",fg_rank," absolute rank",myrank,
 c     & " nfgtasks",nfgtasks
@@ -117,6 +118,9 @@ c        call chainbuild_cart
       endif
 c      write (iout,*) "itime_mat",itime_mat," imatupdate",imatupdate
       if (mod(itime_mat,imatupdate).eq.0) then
+#ifdef TIMING_ENE
+        time01=MPI_Wtime()
+#endif
         call make_SCp_inter_list
 c        write (iout,*) "Finished make_SCp_inter_list"
 c        call flush(iout)
@@ -126,9 +130,12 @@ c        call flush(iout)
         call make_pp_inter_list
 c        write (iout,*) "Finished make_pp_inter_list"
 c        call flush(iout)
-        call make_pp_vdw_inter_list
+c        call make_pp_vdw_inter_list
 c        write (iout,*) "Finished make_pp_vdw_inter_list"
 c        call flush(iout)
+#ifdef TIMING_ENE
+        time_list=time_list+MPI_Wtime()-time01
+#endif
       endif
 c      print *,'Processor',myrank,' calling etotal ipot=',ipot
 c      print *,'Processor',myrank,' nnt=',nnt,' nct=',nct
@@ -151,6 +158,9 @@ C
 C Compute the side-chain and electrostatic interaction energy
 C
 C      print *,ipot
+#ifdef TIMING_ENE
+      time01=MPI_Wtime()
+#endif
       goto (101,102,103,104,105,106) ipot
 C Lennard-Jones potential.
   101 call elj(evdw)
@@ -175,6 +185,9 @@ C
 C Calculate electrostatic (H-bonding) energy of the main chain.
 C
   107 continue
+#ifdef TIMING_ENE
+      time_evdw=time_evdw+MPI_Wtime()-time01
+#endif
 #ifdef DFA
 C     BARTEK for dfa test!
 c      print *,"Processors",MyRank," wdfa",wdfa_dist
@@ -216,6 +229,9 @@ c      print *,"Processor",myrank," computed USCSC"
 #ifdef TIMING
       time_vec=time_vec+MPI_Wtime()-time01
 #endif
+#ifdef TIMING_ENE
+      time01=MPI_Wtime()
+#endif
 C Introduction of shielding effect first for each peptide group
 C the shielding factor is set this factor is describing how each
 C peptide group is shielded by side-chains
@@ -252,6 +268,9 @@ c      print *,"Processor",myrank," left VEC_AND_DERIV"
 c        call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3,
 c     &   eello_turn4)
       endif
+#ifdef TIMING_ENE
+      time_eelec=time_eelec+MPI_Wtime()-time01
+#endif
 c#ifdef TIMING
 c      time_enecalc=time_enecalc+MPI_Wtime()-time00
 c#endif
@@ -260,6 +279,9 @@ C
 C Calculate excluded-volume interaction energy between peptide groups
 C and side chains.
 C
+#ifdef TIMING_ENE
+      time01=MPI_Wtime()
+#endif
       if (ipot.lt.6) then
        if(wscp.gt.0d0) then
         call escp(evdw2,evdw2_14)
@@ -271,6 +293,9 @@ C
 c        write (iout,*) "Soft-sphere SCP potential"
         call escp_soft_sphere(evdw2,evdw2_14)
       endif
+#ifdef TIMING_ENE
+      time_escp=time_escp+MPI_Wtime()-time01
+#endif
 c
 c Calculate the bond-stretching energy
 c
@@ -839,12 +864,12 @@ c      do i=nnt,nres
           gradbufc(k,i)=0.0d0
         enddo
       enddo
-#ifdef DEBUG
-      write (iout,*) "igrad_start",igrad_start," igrad_end",igrad_end
-      write (iout,*) (i," jgrad_start",jgrad_start(i),
-     &                  " jgrad_end  ",jgrad_end(i),
-     &                  i=igrad_start,igrad_end)
-#endif
+c#ifdef DEBUG
+c      write (iout,*) "igrad_start",igrad_start," igrad_end",igrad_end
+c      write (iout,*) (i," jgrad_start",jgrad_start(i),
+c     &                  " jgrad_end  ",jgrad_end(i),
+c     &                  i=igrad_start,igrad_end)
+c#endif
 c
 c Obsolete and inefficient code; we can make the effort O(n) and, therefore,
 c do not parallelize this part.
@@ -1510,6 +1535,7 @@ C
       double precision sscale,sscagrad,sscagradlip,sscalelip
       double precision gg_lipi(3),gg_lipj(3)
       double precision boxshift
+      external boxshift
 c      write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
       gg_lipi=0.0d0
@@ -2039,6 +2065,7 @@ c          do j=istart(i,iint),iend(i,iint)
 c              write(iout,*) "PRZED ZWYKLE", evdwij
               call dyn_ssbond_ene(i,j,evdwij)
 c              write(iout,*) "PO ZWYKLE", evdwij
+c              call flush(iout)
 
               evdw=evdw+evdwij
               if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') 
@@ -2151,7 +2178,7 @@ C I hate to put IF's in the loops, but here don't have another choice!!!!
 cd              write (iout,'(2(a3,i3,2x),17(0pf7.3))')
 cd     &        restyp(itypi),i,restyp(itypj),j,
 cd     &        rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) 
-c                return
+                return
               endif
               sigder=-sig*sigsq
 c---------------------------------------------------------------
@@ -3049,7 +3076,7 @@ c
         write(iout,*)  'b2=',(b2(k,i-2),k=1,2)
 #endif
       enddo
-      mu=0.0d0
+      mu(:,:nres)=0.0d0
 #ifdef PARMAT
       do i=ivec_start+2,ivec_end+2
 #else
@@ -5504,6 +5531,9 @@ C peptide-group centers and side chains and its gradient in virtual-bond and
 C side-chain vectors.
 C
       implicit none
+#ifdef MPI
+      include 'mpif.h'
+#endif
       include 'DIMENSIONS'
       include 'COMMON.GEO'
       include 'COMMON.VAR'
@@ -5515,6 +5545,7 @@ C
       include 'COMMON.IOUNITS'
       include 'COMMON.CONTROL'
       include 'COMMON.SPLITELE'
+      include 'COMMON.TIME1'
       double precision ggg(3)
       integer i,iint,j,k,iteli,itypj,subchap,ikont
       double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1,
@@ -5522,6 +5553,10 @@ C
       double precision evdw2,evdw2_14,evdwij
       double precision sscale,sscagrad
       double precision boxshift
+      external boxshift,to_box
+c#ifdef TIMING_ENE
+c      double precision time01
+c#endif
       evdw2=0.0D0
       evdw2_14=0.0d0
 c        print *,boxxsize,boxysize,boxzsize,'wymiary pudla'
@@ -5533,6 +5568,9 @@ C      do zshift=-1,1
       if (energy_dec) write (iout,*) "escp:",r_cut_int,rlamb
 c      do i=iatscp_s,iatscp_e
       do ikont=g_listscp_start,g_listscp_end
+c#ifdef TIMING_ENE
+c        time01=MPI_Wtime()
+c#endif
         i=newcontlistscpi(ikont)
         j=newcontlistscpj(ikont)
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
@@ -5540,6 +5578,7 @@ c      do i=iatscp_s,iatscp_e
         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))
+!DIR$ INLINE
         call to_box(xi,yi,zi)
 c        do iint=1,nscp_gr(i)
 
@@ -5554,11 +5593,21 @@ C Uncomment following three lines for Ca-p interactions
           xj=c(1,j)
           yj=c(2,j)
           zj=c(3,j)
+!DIR$ INLINE
           call to_box(xj,yj,zj)
+c#ifdef TIMING_ENE
+c       time_escpsetup=time_escpsetup+MPI_Wtime()-time01
+c       time01=MPI_Wtime()
+c#endif
+!DIR$ INLINE
           xj=boxshift(xj-xi,boxxsize)
           yj=boxshift(yj-yi,boxysize)
           zj=boxshift(zj-zi,boxzsize)
 c          print *,xj,yj,zj,'polozenie j'
+c#ifdef TIMING_ENE
+c       time_escpsetup=time_escpsetup+MPI_Wtime()-time01
+c       time01=MPI_Wtime()
+c#endif
           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
 c          print *,rrij
           sss=sscale(1.0d0/(dsqrt(rrij)),r_cut_int)
@@ -5619,6 +5668,9 @@ cgrad          enddo
             gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
             gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
           enddo
+c#ifdef TIMING_ENE
+c          time_escpcalc=time_escpcalc+MPI_Wtime()-time01
+c#endif
 c        endif !endif for sscale cutoff
 c        enddo ! j
 
@@ -5804,7 +5856,8 @@ C 15/02/13 CC dynamic SSbond - additional check
           if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and.
      &        iabs(itype(jjj)).eq.1) then
            call ssbond_ene(iii,jjj,eij)
-           ehpb=ehpb+2*eij
+c           ehpb=ehpb+2*eij
+           ehpb=ehpb+eij
          endif
 cd          write (iout,*) "eij",eij
 cd   &   ' waga=',waga,' fac=',fac
index 34a1bd1..51d5b2d 100644 (file)
@@ -13,7 +13,7 @@ cMS$ATTRIBUTES C ::  proc_proc
 #ifdef MPI
       include "mpif.h"
       double precision weights_(n_ene)
-      double precision time00,time_Bcast,time_BcastW
+      double precision time00
       integer ierror,ierr
 #endif
       include 'COMMON.SETUP'
@@ -29,6 +29,7 @@ cMS$ATTRIBUTES C ::  proc_proc
       include 'COMMON.QRESTR'
       include 'COMMON.MD'
       include 'COMMON.CONTROL'
+      include 'COMMON.TIME1'
       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,
@@ -36,6 +37,9 @@ cMS$ATTRIBUTES C ::  proc_proc
      & esaxs_constr,ehomology_constr,edfator,edfanei,edfabet
       integer i,n_corr,n_corr1
 c      write(iout,'(a,i2)')'Calling etotal_long ipot=',ipot
+#ifdef TIMING_ENE
+      double precision time01
+#endif
       if (modecalc.eq.12.or.modecalc.eq.14) then
 #ifdef MPI
 c        if (fg_rank.eq.0) call int_from_cart1(.false.)
@@ -138,16 +142,25 @@ c      print *,'Processor',myrank,' nnt=',nnt,' nct=',nct
         call MPI_Bcast(itime_mat,1,MPI_INT,king,FG_COMM,IERROR)
       endif
       if (mod(itime_mat,imatupdate).eq.0) then
-        call make_SCp_inter_list
-        call make_SCSC_inter_list
+#ifdef TIMING_ENE
+        time01=MPI_Wtime()
+#endif
+        call make_SCp_inter_list_RESPA
+        call make_SCSC_inter_list_RESPA
         call make_pp_inter_list
-        call make_pp_vdw_inter_list
+        call make_pp_vdw_inter_list_RESPA
+#ifdef TIMING_ENE
+        time_list=time_list+MPI_Wtime()-time01
+#endif
       endif
 #endif     
 cd    print *,'nnt=',nnt,' nct=',nct
 C
 C Compute the side-chain and electrostatic interaction energy
 C
+#ifdef TIMING_ENE
+      time01=MPI_Wtime()
+#endif
       goto (101,102,103,104,105,106) ipot
 C Lennard-Jones potential.
   101 call elj_long(evdw)
@@ -171,8 +184,20 @@ C
 C Calculate electrostatic (H-bonding) energy of the main chain.
 C
   107 continue
+#ifdef TIMING_ENE
+      time_evdw_long=time_evdw_long+MPI_Wtime()-time01
+#endif
+#ifdef TIMING
+      time01=MPI_Wtime() 
+#endif
       call vec_and_deriv
+#ifdef TIMING
+      time_vec=time_vec+MPI_Wtime()-time01
+#endif
 c      write (iout,*) "etotal_long: shield_mode",shield_mode
+#ifdef TIMING_ENE
+      time01=MPI_Wtime()
+#endif
       if (shield_mode.eq.1) then
        call set_shield_fac
       else if  (shield_mode.eq.2) then
@@ -204,10 +229,16 @@ c        write (iout,*) "Soft-spheer ELEC potential"
         call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3,
      &   eello_turn4)
       endif
+#ifdef TIMING_ENE
+      time_eelec_long=time_eelec_long+MPI_Wtime()-time01
+#endif
 C
 C Calculate excluded-volume interaction energy between peptide groups
 C and side chains.
 C
+#ifdef TIMING_ENE
+      time01=MPI_Wtime()
+#endif
       if (ipot.lt.6) then
        if(wscp.gt.0d0) then
         call escp_long(evdw2,evdw2_14)
@@ -218,6 +249,9 @@ C
       else
         call escp_soft_sphere(evdw2,evdw2_14)
       endif
+#ifdef TIMING_ENE
+      time_escp_long=time_escp_long+MPI_Wtime()-time01
+#endif
 #ifdef FOURBODY
 C 
 C 12/1/95 Multi-body terms
@@ -334,12 +368,16 @@ cMS$ATTRIBUTES C ::  proc_proc
       include 'COMMON.CONTROL'
       include 'COMMON.SAXS'
       include 'COMMON.TORCNSTR'
+      include 'COMMON.TIME1'
       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,
      & eliptran,Eafmforce,Etube,
      & esaxs_constr,ehomology_constr,edfator,edfanei,edfabet
       integer i,n_corr,n_corr1
+#ifdef TIMING_ENE
+      double precision time01
+#endif
 c      write(iout,'(a,i2)')'Calling etotal_short ipot=',ipot
 c      call flush(iout)
       if (modecalc.eq.12.or.modecalc.eq.14) then
@@ -432,13 +470,13 @@ C FG slaves receive the WEIGHTS array
           wsaxs=weights(26)
         endif
 c        write (iout,*),"Processor",myrank," BROADCAST weights"
-        call MPI_Bcast(c(1,1),maxres6,MPI_DOUBLE_PRECISION,
+        call MPI_Bcast(c(1,1),6*nres,MPI_DOUBLE_PRECISION,
      &    king,FG_COMM,IERR)
 c        write (iout,*) "Processor",myrank," BROADCAST c"
-        call MPI_Bcast(dc(1,1),maxres6,MPI_DOUBLE_PRECISION,
+        call MPI_Bcast(dc(1,1),6*nres,MPI_DOUBLE_PRECISION,
      &    king,FG_COMM,IERR)
 c        write (iout,*) "Processor",myrank," BROADCAST dc"
-        call MPI_Bcast(dc_norm(1,1),maxres6,MPI_DOUBLE_PRECISION,
+        call MPI_Bcast(dc_norm(1,1),6*nres,MPI_DOUBLE_PRECISION,
      &    king,FG_COMM,IERR)
 c        write (iout,*) "Processor",myrank," BROADCAST dc_norm"
         call MPI_Bcast(theta(1),nres,MPI_DOUBLE_PRECISION,
@@ -470,6 +508,9 @@ c      call int_from_cart1(.false.)
 C
 C Compute the side-chain and electrostatic interaction energy
 C
+#ifdef TIMING_ENE
+      time01=MPI_Wtime()
+#endif
       goto (101,102,103,104,105,106) ipot
 C Lennard-Jones potential.
   101 call elj_short(evdw)
@@ -494,16 +535,31 @@ C
 C Calculate electrostatic (H-bonding) energy of the main chain.
 C
   107 continue
+#ifdef TIMING_ENE
+      time_evdw_short=time_evdw_short+MPI_Wtime()-time01
+#endif
 c
 c Calculate the short-range part of Evdwpp
 c
+#ifdef TIMING_ENE
+      time01=MPI_Wtime()
+#endif
       call evdwpp_short(evdw1)
+#ifdef TIMING_ENE
+      time_eelec_short=time_eelec_short+MPI_Wtime()-time01
+#endif
 c
 c Calculate the short-range part of ESCp
 c
+#ifdef TIMING_ENE
+      time01=MPI_Wtime()
+#endif
       if (ipot.lt.6) then
         call escp_short(evdw2,evdw2_14)
       endif
+#ifdef TIMING_ENE
+      time_escp_short=time_escp_short+MPI_Wtime()-time01
+#endif
 c
 c Calculate the bond-stretching energy
 c
@@ -531,7 +587,13 @@ C energy function
 C
 C Calculate the SC local energy.
 C
+#ifdef TIMING
+      time01=MPI_Wtime() 
+#endif
       call vec_and_deriv
+#ifdef TIMING
+      time_vec=time_vec+MPI_Wtime()-time01
+#endif
       call esc(escloc)
 C
 C Calculate the virtual-bond torsional energy.
index ea009b6..603e6c0 100644 (file)
@@ -13,7 +13,7 @@ C Generate random conformation or chain cut and regrowth.
       logical overlap,back,fail
       integer nstart
       integer i,j,k,it,it1,it2,nit,niter,nsi,maxsi,maxnit
-      double precision gen_theta,gen_phi,dist
+      double precision gen_theta,gen_phi,dist,ran_number
 cd    print *,' CG Processor',me,' maxgen=',maxgen
       maxsi=100
 cd    write (iout,*) 'Gen_Rand_conf: nstart=',nstart
@@ -60,12 +60,16 @@ c       write(iout,*)'theta(3)=',rad2deg*theta(3)
        it1=iabs(itype(i-1))
        it2=iabs(itype(i-2))
        it=iabs(itype(i))
-c       print *,'Gen_Rand_Conf: i=',i,' it=',it,' it1=',it1,' it2=',it2,
-c    &    ' nit=',nit,' niter=',niter,' maxgen=',maxgen
+        if (it.eq.ntyp1 .and. it1.eq.ntyp1) then
+          vbld(i)=ran_number(3.8d0,10.0d0)
+          vbld_inv(i)=1.0d0/vbld(i)
+        endif
+c        write (iout,*) 'Gen_Rand_Conf: i=',i,' it=',it,' it1=',it1,
+c     &    ' it2=',it2,' nit=',nit,' niter=',niter,' maxgen=',maxgen
        phi(i+1)=gen_phi(i+1,it1,it)
        if (back) then
           phi(i)=gen_phi(i+1,it2,it1)
-c         print *,'phi(',i,')=',phi(i)
+c          write(iout,*) 'phi(',i,')=',phi(i)," type",it1,it2,it
          theta(i-1)=gen_theta(it2,phi(i-1),phi(i))
          if (it2.ne.10 .and. it2.ne.ntyp1) then
             nsi=0
@@ -78,7 +82,12 @@ c         print *,'phi(',i,')=',phi(i)
           endif
          call locate_next_res(i-1)
        endif
-       theta(i)=gen_theta(it1,phi(i),phi(i+1))
+        if (it1.ne.ntyp1) then
+         theta(i)=gen_theta(it1,phi(i),phi(i+1))
+        else
+          theta(i)=ran_number(1.326d0,2.548d0)
+        endif
+c        write (iout,*) "i",i," theta",theta(i)
         if (it1.ne.10 .and. it1.ne.ntyp1) then 
         nsi=0
         fail=.true.
@@ -86,10 +95,12 @@ c         print *,'phi(',i,')=',phi(i)
           call gen_side(it1,theta(i),alph(i-1),omeg(i-1),fail)
           nsi=nsi+1
         enddo
+c        write (iout,*) "After forward SC generation:",nsi,maxsi
         if (nsi.gt.maxsi) return1
         endif
        call locate_next_res(i)
        if (overlap(i-1)) then
+c          write (iout,*) "overlap",i-1
          if (nit.lt.maxnit) then
            back=.true.
            nit=nit+1
@@ -107,6 +118,7 @@ c         print *,'phi(',i,')=',phi(i)
            endif
           endif
         else
+c          write (iout,*) "No overlap",i-1
          back=.false.
          nit=0 
          i=i+1
@@ -211,11 +223,55 @@ c--------------------------------------------------------------------------
       double precision function gen_phi(i,it1,it2)
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
+      include 'COMMON.IOUNITS'
       include "COMMON.TORCNSTR"
       include 'COMMON.GEO'
       include 'COMMON.BOUNDS'
-      if (raw_psipred .or. ndih_constr.eq.0) then
+      include 'COMMON.INTERACT'
+      double precision sumprob(3)
+      double precision pinorm
+      external pinorm
+      if (ndih_constr.eq.0) then
         gen_phi=ran_number(-pi,pi) 
+      else if (raw_psipred) then
+        if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
+     &    .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
+          ii=iconstr_dih(i)
+          sumprob(1)=vpsipred(2,ii)
+          sumprob(2)=sumprob(1)+vpsipred(3,ii)
+          sumprob(3)=sumprob(2)+vpsipred(1,ii)
+          aux=ran_number(0.0d0,sumprob(3))
+#ifdef DEBUG
+          write(iout,*)"gen_phi: residue i",i," ii",ii," vpsipred",
+     &      vpsipred(2,ii),vpsipred(3,ii),vpsipred(1,ii)," sumprob",
+     &      sumprob(1),sumprob(2),sumprob(3)
+          write (iout,*) "aux",aux
+#endif
+          if (aux.le.sumprob(1)) then
+#ifdef DEBUG
+            write (iout,*) "hel:",
+     &        (phibound(1,i)-sdihed(1,ndih_constr))*rad2deg,
+     &        (phibound(1,i)+sdihed(1,ndih_constr))*rad2deg
+#endif
+            gen_phi=ran_number(phibound(1,i)-sdihed(1,ndih_constr)
+     &        ,phibound(1,i)+sdihed(1,ndih_constr)) 
+          else if (aux.le.sumprob(2)) then
+#ifdef DEBUG
+            write (iout,*) "ext:",
+     &          (phibound(2,i)-sdihed(2,ndih_constr))*rad2deg,
+     &          (phibound(2,i)+sdihed(2,ndih_constr))*rad2deg
+#endif
+           gen_phi=pinorm(ran_number(phibound(2,i)-sdihed(2,ndih_constr)
+     &        ,phibound(2,i)+sdihed(2,ndih_constr)))
+          else
+#ifdef DEBUG
+            write (iout,*) "coil:",-180.0,180.0
+#endif
+            gen_phi=ran_number(-pi,pi) 
+          endif
+        else
+          gen_phi=ran_number(-pi,pi) 
+        endif
       else
 C 8/13/98 Generate phi using pre-defined boundaries
         gen_phi=ran_number(phibound(1,i),phibound(2,i)) 
@@ -288,7 +344,7 @@ c-------------------------------------------------------------------------
       fail=.false.
       if (the.eq.0.0D0 .or. the.eq.pi) then
 #ifdef MPI
-        write (*,'(a,i4,a,i3,a,1pe14.5)') 
+        write (iout,'(a,i4,a,i3,a,1pe14.5)') 
      & 'CG Processor:',me,' Error in GenSide: it=',it,' theta=',the
 #else
 cd        write (iout,'(a,i3,a,1pe14.5)') 
@@ -778,23 +834,24 @@ c     overlapping residues left, or false otherwise (success)
       include 'COMMON.IOUNITS'
       logical had_overlaps,fail,scfail
       integer ioverlap(maxres),ioverlap_last
+      integer maxit_corr /5000/
 
       had_overlaps=.false.
-      call overlap_sc_list(ioverlap,ioverlap_last)
+      call overlap_sc_list(ioverlap,ioverlap_last,.false.)
       if (ioverlap_last.gt.0) then
         write (iout,*) '#OVERLAPing residues ',ioverlap_last
-        write (iout,'(18i5)') (ioverlap(k),k=1,ioverlap_last)
+        write (iout,'(15i6)') (ioverlap(k),k=1,ioverlap_last)
         had_overlaps=.true.
       endif
 
       maxsi=1000
-      do k=1,1000
+      do k=1,maxit_corr
         if (ioverlap_last.eq.0) exit
 
         do ires=1,ioverlap_last 
           i=ioverlap(ires)
           iti=iabs(itype(i))
-          if (iti.ne.10) then
+          if (iti.ne.10 .and. iti.lt.ntyp1) then
             nsi=0
             fail=.true.
             do while (fail.and.nsi.le.maxsi)
@@ -805,16 +862,15 @@ c     overlapping residues left, or false otherwise (success)
             if(fail) goto 999
           endif
         enddo
-
 c        write (iout,*) "before chaincuild overlap_sc_list: dc0",dc(:,0)
 c        call chainbuild_extconf
 c        write (iout,*) "after chaincuild overlap_sc_list: dc0",dc(:,0)
-        call overlap_sc_list(ioverlap,ioverlap_last)
-        write (iout,*) 'Overlaping residues ',ioverlap_last,
-     &           (ioverlap(j),j=1,ioverlap_last)
+        call overlap_sc_list(ioverlap,ioverlap_last,.false.)
+        write (iout,*)'#Overlaping residues @iter',k,":",ioverlap_last
+        write (iout,*)'Residue list:',(ioverlap(j),j=1,ioverlap_last)
       enddo
 
-      if (k.le.1000.and.ioverlap_last.eq.0) then
+      if (k.le.maxit_corr.and.ioverlap_last.eq.0) then
         scfail=.false.
         if (had_overlaps) then
           write (iout,*) '#OVERLAPing all corrected after ',k,
@@ -823,22 +879,29 @@ c        write (iout,*) "after chaincuild overlap_sc_list: dc0",dc(:,0)
       else
         scfail=.true.
         write (iout,*) '#OVERLAPing NOT all corrected ',ioverlap_last
-        write (iout,'(20i4)') (ioverlap(j),j=1,ioverlap_last)
+        write (iout,'(15i6)') (ioverlap(j),j=1,ioverlap_last)
       endif
 
       return
 
  999  continue
-      write (iout,'(a30,i5,a12,i4)') 
+      write (iout,'(a30,i5,a12,i6)') 
      &               '#OVERLAP FAIL in gen_side after',maxsi,
      &               'iter for RES',i
       scfail=.true.
       return
       end
 
-      subroutine overlap_sc_list(ioverlap,ioverlap_last)
-      implicit real*8 (a-h,o-z)
+      subroutine overlap_sc_list(ioverlap,ioverlap_last,lprn)
+      implicit none
       include 'DIMENSIONS'
+#ifdef MPI
+      include "mpif.h"
+      include "COMMON.SETUP"
+      integer ierror
+      integer ioverlap_last_tab(0:max_fg_procs-1),
+     & ioverlap_all(maxres*max_fg_procs),displs(0:max_fg_procs-1)
+#endif
       include 'COMMON.GEO'
       include 'COMMON.LOCAL'
       include 'COMMON.IOUNITS'
@@ -847,19 +910,47 @@ c        write (iout,*) "after chaincuild overlap_sc_list: dc0",dc(:,0)
       include 'COMMON.FFIELD'
       include 'COMMON.VAR'
       include 'COMMON.CALC'
-      logical fail
+      integer ii,itypi,itypj,itypi1,ind,ikont
+      logical fail,lprn
       integer ioverlap(maxres),ioverlap_last
-      data redfac /0.5D0/
+      double precision redfac /0.5D0/
+      double precision rrij,rij_shift,sig0ij,xi,yi,zi,rcomp,sig
+      double precision dist
 
-      write (iout,*) "overlap_sc_list"
+#ifdef MPI
+      if (nfgtasks.gt.1) then
+        if (fg_rank.eq.0) 
+     &   call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR)
+        call MPI_Bcast(c(1,1),6*nres,MPI_DOUBLE_PRECISION,king,FG_COMM,
+     &    IERROR)
+        call MPI_Bcast(dc(1,0),6*(nres+1),MPI_DOUBLE_PRECISION,king,
+     &   FG_COMM,IERROR)
+        call MPI_Bcast(dc_norm(1,0),6*(nres+1),MPI_DOUBLE_PRECISION,
+     &   king,FG_COMM,IERROR)
+      endif
+#endif
+c      write (iout,*) "overlap_sc_list"
 c      write(iout,*) "iatsc_s",iatsc_s," iatsc_e",iatsc_e
-      write(iout,*) "nnt",nnt," nct",nct
+c      write(iout,*) "nnt",nnt," nct",nct
       ioverlap_last=0
 C Check for SC-SC overlaps and mark residues
 c      print *,'>>overlap_sc nnt=',nnt,' nct=',nct
       ind=0
+#ifdef DEBUG
+      write (iout,*) "FG proecssor",fg_rank," g_listscsc_start",
+     & g_listscsc_start," g_listscsc_end",g_listscsc_end
+      write (*,*) "FG proecssor",fg_rank," g_listscsc_start",
+     & g_listscsc_start," g_listscsc_end",g_listscsc_end
+#endif
 c      do i=iatsc_s,iatsc_e
-      do i=nnt,nct
+      do ikont=g_listscsc_start,g_listscsc_end
+c        write (*,*) "FG processor",fg_rank," loop begins ioverlap_last",
+c     &   ioverlap_last
+        i=newcontlisti(ikont)
+        j=newcontlistj(ikont)
+c      do i=nnt,nct
+c        write (*,*) "FG processor",fg_rank," loop begins ioverlap_last",
+c     &   ioverlap_last,"ikont i j",ikont,i,j
         itypi=iabs(itype(i))
         itypi1=iabs(itype(i+1))
         if (itypi.eq.ntyp1) cycle
@@ -873,7 +964,7 @@ c      do i=iatsc_s,iatsc_e
 c
 c        do iint=1,nint_gr(i)
 c          do j=istart(i,iint),iend(i,iint)
-          do j=i+1,nct
+c          do j=i+1,nct
             ind=ind+1
             itypj=iabs(itype(j))
             if (itypj.eq.ntyp1) cycle
@@ -919,10 +1010,16 @@ c     &       " sig",sig," sig0ij",sig0ij
 c            if ( rij_shift.le.0.0D0 ) then
             if ( rij_shift/sig0ij.le.0.1D0 ) then
 c              write (iout,*) "overlap",i,j
-              write (iout,'(a,i5,a,i5,a,f10.5,a,3f10.5)')
-     &         'overlap SC-SC: i=',i,' j=',j,
-     &         ' dist=',dist(nres+i,nres+j),' rcomp=',
-     &         rcomp,1.0/rij,rij_shift
+              if (lprn) then
+                write (iout,'(a,i5,a,i5,a,f10.5,a,3f10.5)')
+     &           'overlap SC-SC: i=',i,' j=',j,
+     &           ' dist=',dist(nres+i,nres+j),' rcomp=',
+     &           rcomp,1.0/rij,rij_shift
+                write (*,'(a,i2,a,i5,a,i5,a,f10.5,a,3f10.5)')
+     &           'FG processor',fg_rank,' overlap SC-SC: i=',i,' j=',j,
+     &           ' dist=',dist(nres+i,nres+j),' rcomp=',
+     &           rcomp,1.0/rij,rij_shift
+              endif
               ioverlap_last=ioverlap_last+1
               ioverlap(ioverlap_last)=i         
               do k=1,ioverlap_last-1
@@ -933,9 +1030,73 @@ c              write (iout,*) "overlap",i,j
               do k=1,ioverlap_last-1
                 if (ioverlap(k).eq.j) ioverlap_last=ioverlap_last-1
               enddo 
+c              write(*,*) "FG processor",fg_rank,i,j," ioverlap_last",
+c     &        ioverlap_last," ioverlap",(ioverlap(k),k=1,ioverlap_last)
             endif
-          enddo
+c          enddo
 c        enddo
       enddo
+#ifdef MPI
+#ifdef DEBUG
+      write (iout,*) "FG Processor",fg_rank," ioverlap_last",
+     & ioverlap_last," ioverlap",(ioverlap(i),i=1,ioverlap_last)
+      write (*,*) "FG Processor",fg_rank," ioverlap_last",ioverlap_last,
+     & " ioverlap",(ioverlap(i),i=1,ioverlap_last)
+      call flush(iout)
+#endif
+      if (nfgtasks.eq.1) return
+#ifdef DEBUG
+      write (iout,*) "Before MPI_Gather"
+      call flush(iout)
+#endif
+      call MPI_Gather(ioverlap_last,1,MPI_INTEGER,ioverlap_last_tab,
+     & 1,MPI_INTEGER,king,FG_COMM,IERROR)
+#ifdef DEBUG
+      write (iout,*) "After MPI_Gather"
+      call flush(iout)
+#endif
+#ifdef DEBUG
+      if (myrank.eq.king) 
+     & write (iout,*) "FG Processor",fg_rank,"ioverlap_last_tab",
+     & (ioverlap_last_tab(i),i=0,nfgtasks-1)
+      call flush(iout)
+#endif
+      displs(0)=0
+      do i=1,nfgtasks-1
+        displs(i)=displs(i-1)+ioverlap_last_tab(i-1)
+      enddo
+      call MPI_Gatherv(ioverlap,ioverlap_last,MPI_INTEGER,
+     & ioverlap_all,ioverlap_last_tab,displs,MPI_INTEGER,king,
+     & FG_COMM,IERROR)
+#ifdef DEBUG
+      write (iout,*) "After Gatherv"
+      call flush(iout)
+#endif
+      if (fg_rank.gt.0) return
+      ioverlap_last=0
+      do i=0,nfgtasks-1
+        ioverlap_last=ioverlap_last+ioverlap_last_tab(i)
+      enddo
+#ifdef DEBUG
+      write (iout,*) "ioverlap_last",ioverlap_last," ioverlap_last",
+     & (ioverlap_all(i),i=1,ioverlap_last)
+      call flush(iout)
+#endif
+      ii=0
+      do i=1,ioverlap_last
+        ioverlap(ii+1)=ioverlap_all(i)
+        do j=ii,1,-1
+          if (ioverlap(ii+1).eq.ioverlap(j)) goto 11
+        enddo
+        ii=ii+1
+   11   continue 
+      enddo
+      ioverlap_last=ii
+#ifdef DEBUG
+      write (iout,*) "After summing: ioverlap_last",ioverlap_last,
+     & " ioverlap",(ioverlap(i),i=1,ioverlap_last)
+      call flush(iout)
+#endif
+#endif
       return
       end
diff --git a/source/unres/src-HCD-5D/gen_rand_conf_mchain.F b/source/unres/src-HCD-5D/gen_rand_conf_mchain.F
new file mode 100644 (file)
index 0000000..4614e87
--- /dev/null
@@ -0,0 +1,424 @@
+      subroutine gen_rand_conf_mchain(nstart0,*)
+C Generate random conformation or chain cut and regrowth.
+      implicit none
+      include 'DIMENSIONS'
+      include 'COMMON.CHAIN'
+      include 'COMMON.LOCAL'
+      include 'COMMON.VAR'
+      include 'COMMON.INTERACT'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.MCM'
+      include 'COMMON.GEO'
+      include 'COMMON.CONTROL'
+      logical overlap_mchain,back,fail
+      integer nstart,nstart0
+      integer i,ii,iii,j,k,it,it1,it2,nit,niter,nsi,maxsi,maxnit
+      integer igr,iequi,ichain,nnres
+      double precision aux
+      double precision gen_theta,gen_phi,dist,ran_number,scalar
+c      write (iout,*)  'gen_rand_conf_mchain: maxgen=',maxgen
+      nstart=nstart0
+      maxsi=100
+c      write (iout,*) 'Gen_Rand_conf_mchain: nstart=',nstart,
+c     &   " nchain_group",nchain_group
+
+      DO IGR=1,NCHAIN_GROUP
+
+      DO IEQUI=1,NEQUIV(IGR)
+
+      ichain=iequiv(iequi,igr)
+
+      i=chain_border1(1,ichain)+nstart-1
+      if (nstart.eq.1) then
+        do j=1,3
+          c(j,i)=ran_number(-15.0d0,15.0d0)
+          dc(j,i-1)=c(j,i)
+        enddo
+      endif
+      if (nstart.le.2) then
+
+      do j=1,3
+        dc_norm(j,i)=ran_number(-1.0d0,1.0d0)
+      enddo
+      aux=dsqrt(scalar(dc_norm(1,i),dc_norm(1,i)))
+      do j=1,3
+        dc_norm(j,i)=dc_norm(j,i)/aux
+      enddo
+      if (itype(i).eq.ntyp1) then
+        do j=1,3
+          dc(j,i)=1.9d0*dc_norm(j,i)
+        enddo
+      else 
+        do j=1,3
+          dc(j,i)=3.8d0*dc_norm(j,i)
+        enddo
+      endif
+      do j=1,3
+        c(j,i+1)=c(j,i)+dc(j,i)
+      enddo
+      endif
+c      if (nstart.lt.5) then
+      if (nstart.le.2) then
+      it1=iabs(itype(i+2))
+      phi(i+3)=gen_phi(i+3,iabs(itype(i+1)),iabs(itype(i+2)))
+c       write(iout,*)'phi(4)=',rad2deg*phi(4)
+      theta(i+2)=gen_theta(iabs(itype(i+2)),pi,phi(i+3))
+      if (it1.ne.10) then
+        nsi=0
+        fail=.true.
+        do while (fail.and.nsi.le.maxsi)
+          call gen_side(it1,theta(i+2),alph(i+1),omeg(i+1),fail)
+          nsi=nsi+1
+        enddo
+        if (nsi.gt.maxsi) then
+          write (iout,'(a,i7,a,i7,a,i7,a,i7)') 
+     &    'Problem in SC rotamer generation, residue',ii,
+     &    ' chain',ichain,' chain group',igr,'. Increase MAXSI.'
+          return1
+        endif
+      endif ! it1.ne.10
+      call orig_frame_chain(i)
+      else
+c      call orig_frame_chain(i-1)
+c      write (iout,*) "calling refsys",i
+      call refsys(i,i-1,i-2,prod(1,1,i-1),
+     &  prod(1,2,i-1),prod(1,3,i-1),fail)
+c      write (iout,*) "after refsys",i
+#ifdef DEBUG
+      write (iout,*) "dc_norm(:",i-1,") and prod"
+      do j=1,3
+        write (iout,*) j,dc_norm(j,i-1),(prod(j,k,i-1),k=1,3)
+      enddo
+#endif
+      endif
+
+      ENDDO
+
+      nstart=nstart+1
+      ii=nstart
+
+      maxnit=5000
+
+      nit=0
+      niter=0
+      back=.false.
+      nnres=chain_border1(2,iequiv(1,igr))-
+     &  chain_border1(1,iequiv(1,igr))+1
+#ifdef DEBUG
+      write (iout,*) "chain group",igr," chains",
+     &    (iequiv(j,igr),j=1,nequiv(igr))
+      write (iout,*) "ii",ii," nnres",nnres
+#endif
+      do while (ii.le. nnres .and. niter.lt.maxgen)
+
+        ichain=iequiv(1,igr)
+        i=ii-1+chain_border1(1,ichain)
+#ifdef DEBUG
+        write (iout,*) "ii",ii," nnres",nnres," ichain",ichain," i",i,
+     &    "niter",niter," back",back," nstart",nstart
+#endif
+        if (ii.lt.nstart) then
+c          if(iprint.gt.1) then
+          write (iout,'(/80(1h*)/2a/80(1h*))') 
+     &          'Generation procedure went down to ',
+     &          'chain beginning. Cannot continue...'
+          write (*,'(/80(1h*)/2a/80(1h*))') 
+     &          'Generation procedure went down to ',
+     &          'chain beginning. Cannot continue...'
+c          endif
+          return1
+        endif
+        it1=iabs(itype(i-1))
+        it2=iabs(itype(i-2))
+        it=iabs(itype(i))
+        if (it.eq.ntyp1 .and. it1.eq.ntyp1) then
+          vbld(i)=ran_number(3.8d0,10.0d0)
+          vbld_inv(i)=1.0d0/vbld(i)
+        endif
+#ifdef DEBUG
+        write (iout,*) 'Gen_Rand_Conf: i=',i,' it=',it,' it1=',it1,
+     &    ' it2=',it2,' nit=',nit,' niter=',niter,' maxgen=',maxgen
+#endif
+        phi(i)=gen_phi(i,it2,it1)
+        phi(i+1)=gen_phi(i+1,it1,it)
+#ifdef DEBUG
+        write (iout,*) "phi",i,phi(i)," phi",i+1,phi(i+1)
+#endif
+        do iequi=2,nequiv(igr)
+          iii=ii+chain_border1(1,iequiv(iequi,igr))
+          phi(iii)=phi(i+1)
+        enddo
+#ifdef CHUJ
+        if (back) then
+          phi(i)=gen_phi(i+1,it2,it1)
+#ifdef DEBUG
+          write(iout,*) 'phi(',i,')=',phi(i)," type",it1,it2,it
+#endif
+          theta(i-1)=gen_theta(it2,phi(i-1),phi(i))
+          if (theta(i-1).gt.2.68780478d0) theta(i-1)=2.68780478d0
+          do iequi=2,nequiv(igr)
+            iii=ii-1+chain_border1(1,iequiv(iequi,igr))
+            phi(iii)=phi(i+1)
+            theta(iii-1)=theta(i-1)
+          enddo
+          if (it2.ne.10 .and. it2.ne.ntyp1) then
+            nsi=0
+            fail=.true.
+            do while (fail.and.nsi.le.maxsi)
+              call gen_side(it2,theta(i-1),alph(i-2),omeg(i-2),fail)
+              nsi=nsi+1
+            enddo
+            do iequi=2,nequiv(igr)
+              iii=ii-3+chain_border1(1,iequiv(iequi,igr))
+              alph(iii)=alph(i-2)
+              omeg(iii)=omeg(i-2)
+            enddo
+#ifdef DEBUG
+            write (iout,*) "alpha",alph(i)," omeg",omeg(i)," fail",fail
+#endif
+            if (nsi.gt.maxsi) then
+              write (iout,'(a,i7,a,i7,a,i7,a,i7)') 
+     &        'Problem in SC rotamer generation, residue',ii,
+     &        ' chain',ichain,' chain group',igr,'. Increase MAXSI.'
+              return1
+            endif
+
+          endif
+c          call locate_next_res(i-1)
+        endif
+#endif
+        if (it1.ne.ntyp1) then
+          theta(i)=gen_theta(it1,phi(i),phi(i+1))
+          if (theta(i).gt.2.68780478d0) theta(i)=2.68780478d0
+        else
+          theta(i)=ran_number(1.326d0,2.548d0)
+        endif
+#ifdef DEBUG
+        write (iout,*) "ii",ii," i",i," it1",it1,
+     &   " theta",theta(i)," phi",
+     &    phi(i),phi(i+1)
+#endif
+        if (it1.ne.10 .and. it1.ne.ntyp1) then 
+        nsi=0
+        fail=.true.
+        do while (fail.and.nsi.le.maxsi)
+c          theta(i)=gen_theta(it1,phi(i),phi(i+1))
+          call gen_side(it1,theta(i),alph(i-1),omeg(i-1),fail)
+          nsi=nsi+1
+        enddo
+#ifdef DEBUG
+        write (iout,*) "alpha",alph(i)," omeg",omeg(i)," fail",fail
+#endif
+c        write (iout,*) "After forward SC generation:",nsi,maxsi
+        if (nsi.gt.maxsi) then
+          write (iout,'(a,i7,a,i7,a,i7,a,i7)') 
+     &    'Problem in SC rotamer generation, residue',ii,
+     &    ' chain',ichain,' chain group',igr,'. Increase MAXSI.'
+          return1
+        endif
+
+        endif
+
+        do iequi=2,nequiv(igr)
+          iii=ii-1+chain_border1(1,iequiv(iequi,igr))
+          theta(iii)=theta(i)
+          phi(iii)=phi(i)
+          alph(iii-1)=alph(i-1)
+          omeg(iii-1)=omeg(i-1)
+        enddo
+
+        DO IEQUI=1,NEQUIV(IGR)
+
+        ichain=iequiv(iequi,igr)
+        i=ii-1+chain_border1(1,ichain)
+#ifdef CHUJ
+        if (back) call locate_next_res(i-1)
+#endif
+        call locate_next_res(i)
+#ifdef DEBUG
+        write (iout,*) "i",i," ii",ii," ichain",ichain
+        write (iout,*) theta(i)*rad2deg,phi(i)*rad2deg,
+     &    alph(i-1)*rad2deg,omeg(i-1)*rad2deg
+        write (iout,*) (c(j,i),j=1,3),(c(j,i+nres-1),j=1,3)
+#endif
+        if (overlap_mchain(i-1,ii-1,ichain,igr)) then
+#ifdef DEBUG
+          write (iout,*) "***********overlap",i-1," nit",nit
+#endif
+          if (nit.lt.maxnit) then
+            back=.true.
+            nit=nit+1
+            exit
+          else
+#ifdef DEBUG
+            write (iout,*) "***********overlap maxnit exceeded",nit
+#endif
+            nit=0
+            if (ii.gt.3) then
+              back=.true.
+              ii=ii-1
+c              write (iout,*) "ii",ii
+              exit
+            else
+              write (iout,'(a,i7,a,i7,a,i7,a,i7)') 
+     &        'Cannot generate non-overlaping conformation, residue',ii,
+     &        ' chain',ichain,' chain group',igr,'. Increase MAXNIT.'
+              return1
+            endif
+          endif
+        else
+c          write (iout,*) "No overlap",i-1
+          back=.false.
+c          nit=0 
+c          i=i+1
+        endif
+
+        ENDDO
+
+        if (.not.back) then
+#ifdef DEBUG
+          write (iout,*) "++++++++++Successful generation",igr,ichain,ii
+#endif
+          ii=ii+1
+          nit=0
+        endif
+        back=.false.
+c        write (iout,*) "ii",ii
+        niter=niter+1
+
+      enddo
+      if (niter.ge.maxgen) then
+        write (iout,'(a,2i7,a,i7,a,i7,a,i7)') 
+     & 'Too many trials in conformation generation',niter,maxgen,
+     & ' chain group',igr,' chain',ichain,' residue',ii
+        return1
+      endif
+
+      ENDDO
+
+      do i=2,nres
+        if (itype(i).eq.ntyp1) then
+          do j=1,3
+            dc(j,i-1)=c(j,i)-c(j,i-1)
+          enddo
+        endif
+      enddo
+
+      return
+      end
+c-------------------------------------------------------------------------
+      logical function overlap_mchain(i,ii,ichain,igr)
+      implicit none
+      include 'DIMENSIONS'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CHAIN'
+      include 'COMMON.INTERACT'
+      include 'COMMON.FFIELD'
+      double precision redfac /0.5D0/,redfacp /0.8d0/,redfacscp /0.8d0/
+      integer i,ii,ichain,igr,j,jj,jchain,jgr,k,iti,itj,iteli,itelj
+      double precision rcomp
+      double precision dist
+      logical lprn /.false./
+      overlap_mchain=.false.
+      iti=iabs(itype(i))
+      if (iti.gt.ntyp) return
+C Check for SC-SC overlaps.
+cd    print *,'nnt=',nnt,' nct=',nct
+#ifdef DEBUG
+      write (iout,*) "overlap_mchain i",i," ii",ii," ichain",ichain,
+     &   " igr",igr," iti",iti
+#endif
+      do j=nnt,nct
+        if (itype(j).eq.ntyp1) cycle
+        jchain=ireschain(j)
+c        write (iout,*) "overlap_mchain j",j," jj",jj," jchain",jchain,
+c     &    " itype",itype(j)
+        jgr=mapchain(jchain)
+        jj=j-chain_border1(1,jchain)+1
+        itj=iabs(itype(j))
+c        write (iout,*) "jgr",jgr
+        if(igr.eq.jgr.and.jj.gt.ii .or. ichain.eq.jchain .and. j.gt.i-1
+     &    .or. jgr.gt.igr .and. jj.gt.nran_start) cycle
+        if (j.lt.i-1 .or. ipot.ne.4) then
+          rcomp=sigmaii(iti,itj)
+        else 
+          rcomp=sigma(iti,itj)
+        endif
+cd      print *,'j=',j
+        if (dist(nres+i,nres+j).lt.redfac*rcomp) then
+          overlap_mchain=.true.
+          if (lprn) write(iout,*)'overlap_mchain, SC-SC: i=',i,' j=',j,
+     &     ' ichain',ichain,' jchain',jchain,
+     &     ' dist=',dist(nres+i,nres+j),' rcomp=',
+     &     rcomp*redfac
+        return
+        endif
+      enddo
+#ifdef CHUJ
+C Check for overlaps between the added peptide group and the preceding
+C SCs.
+      iteli=itel(i)
+      if (iteli.gt.0) then
+      do j=1,3
+        c(j,maxres2+1)=0.5D0*(c(j,i)+c(j,i+1))
+      enddo
+      do j=nnt,nct
+        itj=iabs(itype(j))
+        if (itj.eq.ntyp1) cycle
+        jchain=ireschain(j)
+        jgr=mapchain(jchain)
+        jj=j-chain_border1(1,jchain)+1
+        if(igr.eq.jgr.and.jj.gt.ii .or. ichain.eq.jchain .and. j.gt.i-2
+     &    .or. jgr.gt.igr .and. jj.gt.nran_start) cycle
+cd      print *,'overlap_mchain, p-Sc: i=',i,' j=',j,
+cd   &         ' dist=',dist(nres+j,maxres2+1)
+        if (dist(j,maxres2+1).lt.4.0D0*redfacscp) then
+          if (lprn) write (iout,*) 'overlap_mchain, p-SC: i=',i,' j=',j,
+     &     ' ichain',ichain,' jchain',jchain,
+     &     ' dist=',dist(nres+j,maxres2+1),' rcomp=',
+     &    4.0d0*redfac 
+          overlap_mchain=.true.
+          return
+        endif
+      enddo
+      endif
+C Check for p-p overlaps
+#endif
+      iteli=itel(i)
+      if (iteli.eq.0) return
+      do j=1,3
+        c(j,maxres2+1)=0.5D0*(c(j,i)+c(j,i+1))
+      enddo
+c      do j=nnt,i-2
+      do j=nnt,nct
+        itelj=itel(j)
+        if (itelj.eq.0) cycle
+        do k=1,3
+          c(k,maxres2+2)=0.5D0*(c(k,j)+c(k,j+1))
+        enddo
+        jchain=ireschain(j)
+        jgr=mapchain(jchain)
+c        if(iteli.ne.0.and.itelj.ne.0)then
+c        write (iout,*) i,j,dist(maxres2+1,maxres2+2),rpp(iteli,itelj)
+        jj=j-chain_border1(1,jchain)+1
+c        write (iout,*) "jgr",jgr
+        if(igr.eq.jgr.and.jj.gt.ii .or. ichain.eq.jchain .and. j.gt.i-2
+     &    .or. jgr.gt.igr .and. jj.gt.nran_start) cycle
+#ifdef DEBUG
+          write (iout,*)'overlap_mchain, p-p: i=',i,' j=',j,' igr',igr,
+     &   ' jgr',jgr,' ichain',ichain,' jchain',jchain,
+     &   ' dist=',dist(maxres2+1,maxres2+2)
+#endif
+        if (dist(maxres2+1,maxres2+2).lt.rpp(iteli,itelj)*redfacp) then
+          if (lprn) write (iout,*) 'overlap_mchain, p-p: i=',i,' j=',j,
+     &     ' ichain',ichain,' jchain',jchain,
+     &     ' dist=',dist(maxres2+1,maxres2+2),' rcomp=',
+     &    rpp(iteli,itelj)*redfacp
+          overlap_mchain=.true.
+          return
+        endif
+c        endif
+      enddo
+      return
+      end
+
index d1a3a87..3dcde10 100644 (file)
 #endif
       character*50 tytul
       integer iunit
-      character*1 chainid(10) /'A','B','C','D','E','F','G','H','I','J'/
+      character*1 chainid(52) /'A','B','C','D','E','F','G','H','I','J',
+     & 'K','L','M','N','O','P','Q','R','S','T','U','V','W','X','Y','Z',
+     & 'a','b','c','d','e','f','g','h','i','j','k','l','m','n','o','p',
+     & 'q','r','s','t','u','v','w','x','y','z'/
       integer ica(maxres)
       integer i,j,k,iti,itj,itk,itl,iatom,ichain,ires
       double precision etot
@@ -90,8 +93,8 @@ cmodel      write (iunit,'(a5,i6)') 'MODEL',1
         do i=1,nss
          if (dyn_ss) then
           write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') 
-     &         'SSBOND',i,'CYS',idssb(i)-nnt+1,
-     &                    'CYS',jdssb(i)-nnt+1
+     &         'SSBOND',i,'CYS',iss(idssb(i))-nnt+1,
+     &                    'CYS',iss(jdssb(i))-nnt+1
          else
           write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') 
      &         'SSBOND',i,'CYS',ihpb(i)-nnt+1-nres,
@@ -107,9 +110,10 @@ cmodel      write (iunit,'(a5,i6)') 'MODEL',1
         iti=itype(i)
         if ((iti.eq.ntyp1).and.((itype(i+1)).eq.ntyp1)) then
           ichain=ichain+1
+          if (ichain.gt.52) ichain=1
           ires=0
           write (iunit,'(a)') 'TER'
-        else
+        else if (iti.ne.ntyp1) then
         ires=ires+1
         iatom=iatom+1
         ica(i)=iatom
@@ -149,7 +153,7 @@ cmodel      write (iunit,'(a5,i6)') 'MODEL',1
       write (iunit,'(a6)') 'ENDMDL'     
   10  FORMAT ('ATOM',I7,'  CA  ',A3,1X,A1,I4,4X,3F8.3,f15.3)
   20  FORMAT ('ATOM',I7,'  CB  ',A3,1X,A1,I4,4X,3F8.3,f15.3)
-  30  FORMAT ('CONECT',8I5)
+  30  FORMAT ('CONECT',8I7)
       return
       end
 c------------------------------------------------------------------------------
@@ -541,7 +545,7 @@ C          print *,'A CHUJ',potEcomp(23)
           if(itime.eq.0) then
            write(format,'(a1,a4,a1,a4,a10)') "(",format1,",",format2,
      &                                                     ",31a12)"
-           write (istat,format) "#","",
+           write (istat,format) "#"," ",
      &      (ename(print_order(i)),i=1,nprint_ene)
           endif
           write(format,'(a1,a4,a1,a4,a10)') "(",format1,",",format2,
index 67275ed..af7978b 100644 (file)
@@ -210,6 +210,7 @@ C-------------------------------------------------------------------------
       include 'COMMON.IOUNITS'
       include 'COMMON.TIME1'
       integer i,j,kk
+      double precision time00,time01
 c
 c This subrouting calculates total Cartesian coordinate gradient. 
 c The subroutine chainbuild_cart and energy MUST be called beforehand.
index 710f907..6a297b8 100644 (file)
@@ -351,7 +351,7 @@ c-------------------------------------------------------------------------
      1   "WSC   ","WSCP  ","WELEC","WCORR","WCORR5","WCORR6","WEL_LOC",
 !          8        9       10      11      12       13       14
      8   "WTURN3","WTURN4","WTURN6","WANG","WSCLOC","WTOR  ","WTORD  ",
-!         15       16       17      18      19       20       21
+!         15       16       17      18      19       20       21
      5   "WSTRAIN","WVDWPP","WBOND","SCAL14","WDIHC","WUMB","WSCCOR",
 !         22       23       24      25      26       27       28
      2   "WLT","WAFM","WTHETCNSR","WTUBE","WSAXS","WHOMO","WDFAD",
@@ -453,8 +453,13 @@ c---------------------------------------------------------------------------
       integer iturn3_start_all,iturn3_end_all,iturn4_start_all,
      & iturn4_end_all,iatel_s_all,
      & iatel_e_all,ielstart_all,ielend_all,ntask_cont_from_all,
-     & itask_cont_from_all,ntask_cont_to_all,itask_cont_to_all,
-     & n_sc_int_tot,my_sc_inds,my_sc_inde,ind_sctint,ind_scint_old
+     & itask_cont_from_all,ntask_cont_to_all,itask_cont_to_all
+      integer*8 n_sc_int_tot,my_sc_inds,my_sc_inde,ind_scint,
+     & ind_scint_old,nele_int_tot,ind_eleint,my_ele_inds,my_ele_inde,
+     & ind_eleint_old,nele_int_tot_vdw,my_ele_inds_vdw,my_ele_inde_vdw,
+     & ind_eleint_vdw,ind_eleint_vdw_old,nscp_int_tot,my_scp_inds,
+     & my_scp_inde,ind_scpint,ind_scpint_old,nsumgrad,nlen,ngrad_start,
+     & ngrad_end
       common /przechowalnia/ iturn3_start_all(0:max_fg_procs),
      & iturn3_end_all(0:max_fg_procs),iturn4_start_all(0:max_fg_procs),
      & iturn4_end_all(0:max_fg_procs),iatel_s_all(0:max_fg_procs),
@@ -466,12 +471,8 @@ c---------------------------------------------------------------------------
      & itask_cont_to_all(0:max_fg_procs-1,0:max_fg_procs-1)
       integer FG_GROUP,CONT_FROM_GROUP,CONT_TO_GROUP
       logical scheck,lprint,flag
-      integer i,j,k,ii,jj,iint,npept,nele_int_tot,ind_eleint,ind_scint,
-     & my_ele_inds,my_ele_inde,ind_eleint_old,nele_int_tot_vdw,
-     & my_ele_inds_vdw,my_ele_inde_vdw,ind_eleint_vdw,ijunk,
-     & ind_eleint_vdw_old,nscp_int_tot,my_scp_inds,my_scp_inde,
-     & ind_scpint,ind_scpint_old,nsumgrad,nlen,ngrad_start,ngrad_end,
-     & iaux,ind_typ,ncheck_from,ncheck_to,ichunk
+      integer i,j,k,ii,jj,iint,npept,
+     & ijunk,iaux,ind_typ,ncheck_from,ncheck_to,ichunk
 #ifdef MPI
       integer my_sc_int(0:max_fg_Procs-1),my_sc_intt(0:max_fg_Procs),
      & my_ele_int(0:max_fg_Procs-1),my_ele_intt(0:max_fg_Procs)
@@ -486,14 +487,14 @@ C... to deal with by current processor.
       lprint=energy_dec
       if (lprint)
      &write (iout,*) 'INIT_INT_TABLE nres=',nres,' nnt=',nnt,' nct=',nct
-      n_sc_int_tot=(nct-nnt+1)*(nct-nnt)/2-nss
-      call int_bounds(n_sc_int_tot,my_sc_inds,my_sc_inde)
+      n_sc_int_tot=int(nct-nnt+1,8)*int(nct-nnt,8)/2-nss
+      call int_bounds8(n_sc_int_tot,my_sc_inds,my_sc_inde)
       if (lprint)
      &  write (iout,*) 'Processor',fg_rank,' CG group',kolor,
      &  ' absolute rank',MyRank,
      &  ' n_sc_int_tot',n_sc_int_tot,' my_sc_inds=',my_sc_inds,
      &  ' my_sc_inde',my_sc_inde
-      ind_sctint=0
+      ind_scint=0
       iatsc_s=0
       iatsc_e=0
 #endif
@@ -530,7 +531,7 @@ cd      write (iout,*) 'i=',i,' scheck=',scheck,' jj=',jj
           if (jj.eq.i+1) then
 #ifdef MPI
 c            write (iout,*) 'jj=i+1'
-            call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
+            call int_partition8(ind_scint,my_sc_inds,my_sc_inde,i,
      & iatsc_s,iatsc_e,i+2,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
 #else
             nint_gr(i)=1
@@ -540,7 +541,7 @@ c            write (iout,*) 'jj=i+1'
           else if (jj.eq.nct) then
 #ifdef MPI
 c            write (iout,*) 'jj=nct'
-            call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
+            call int_partition8(ind_scint,my_sc_inds,my_sc_inde,i,
      &  iatsc_s,iatsc_e,i+1,nct-1,nint_gr(i),istart(i,1),iend(i,1),*12)
 #else
             nint_gr(i)=1
@@ -549,10 +550,10 @@ c            write (iout,*) 'jj=nct'
 #endif
           else
 #ifdef MPI
-            call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
+            call int_partition8(ind_scint,my_sc_inds,my_sc_inde,i,
      & iatsc_s,iatsc_e,i+1,jj-1,nint_gr(i),istart(i,1),iend(i,1),*12)
             ii=nint_gr(i)+1
-            call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
+            call int_partition8(ind_scint,my_sc_inds,my_sc_inde,i,
      & iatsc_s,iatsc_e,jj+1,nct,nint_gr(i),istart(i,ii),iend(i,ii),*12)
 #else
             nint_gr(i)=2
@@ -564,7 +565,7 @@ c            write (iout,*) 'jj=nct'
           endif
         else
 #ifdef MPI
-          call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
+          call int_partition8(ind_scint,my_sc_inds,my_sc_inde,i,
      &    iatsc_s,iatsc_e,i+1,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
 #else
           nint_gr(i)=1
@@ -588,9 +589,10 @@ c            write (iout,*) 'jj=nct'
      &   ' absolute rank',myrank,' iatsc_s=',iatsc_s,' iatsc_e=',iatsc_e
 #endif
       if (lprint) then
+      write (iout,*) 'iatsc_s=',iatsc_s,' iatsc_e=',iatsc_e
       write (iout,'(a)') 'Interaction array:'
       do i=iatsc_s,iatsc_e
-        write (iout,'(i3,2(2x,2i3))') 
+        write (iout,'(i7,2(2x,2i7))') 
      & i,(istart(i,iint),iend(i,iint),iint=1,nint_gr(i))
       enddo
       endif
@@ -598,8 +600,8 @@ c            write (iout,*) 'jj=nct'
 #ifdef MPI
 C Now partition the electrostatic-interaction array
       npept=nct-nnt
-      nele_int_tot=(npept-ispp)*(npept-ispp+1)/2
-      call int_bounds(nele_int_tot,my_ele_inds,my_ele_inde)
+      nele_int_tot=int(npept-ispp,8)*int(npept-ispp+1,8)/2
+      call int_bounds8(nele_int_tot,my_ele_inds,my_ele_inde)
       if (lprint)
      & write (*,*) 'Processor',fg_rank,' CG group',kolor,
      &  ' absolute rank',MyRank,
@@ -611,14 +613,14 @@ C Now partition the electrostatic-interaction array
       ind_eleint_old=0
       do i=nnt,nct-3
         ijunk=0
-        call int_partition(ind_eleint,my_ele_inds,my_ele_inde,i,
+        call int_partition8(ind_eleint,my_ele_inds,my_ele_inde,i,
      &    iatel_s,iatel_e,i+ispp,nct-1,ijunk,ielstart(i),ielend(i),*13)
       enddo ! i 
    13 continue
       if (iatel_s.eq.0) iatel_s=1
-      nele_int_tot_vdw=(npept-2)*(npept-2+1)/2
+      nele_int_tot_vdw=int(npept-2,8)*int(npept-2+1,8)/2
 c      write (iout,*) "nele_int_tot_vdw",nele_int_tot_vdw
-      call int_bounds(nele_int_tot_vdw,my_ele_inds_vdw,my_ele_inde_vdw)
+      call int_bounds8(nele_int_tot_vdw,my_ele_inds_vdw,my_ele_inde_vdw)
 c      write (iout,*) "my_ele_inds_vdw",my_ele_inds_vdw,
 c     & " my_ele_inde_vdw",my_ele_inde_vdw
       ind_eleint_vdw=0
@@ -627,7 +629,7 @@ c     & " my_ele_inde_vdw",my_ele_inde_vdw
       iatel_e_vdw=0
       do i=nnt,nct-3
         ijunk=0
-        call int_partition(ind_eleint_vdw,my_ele_inds_vdw,
+        call int_partition8(ind_eleint_vdw,my_ele_inds_vdw,
      &    my_ele_inde_vdw,i,
      &    iatel_s_vdw,iatel_e_vdw,i+2,nct-1,ijunk,ielstart_vdw(i),
      &    ielend_vdw(i),*15)
@@ -655,15 +657,15 @@ c     &   " ielend_vdw",ielend_vdw(i)
      &  ' absolute rank',MyRank
         write (iout,*) 'Electrostatic interaction array:'
         do i=iatel_s,iatel_e
-          write (iout,'(i3,2(2x,2i3))') i,ielstart(i),ielend(i)
+          write (iout,'(i7,2(2x,2i7))') i,ielstart(i),ielend(i)
         enddo
       endif ! lprint
 c     iscp=3
       iscp=2
 C Partition the SC-p interaction array
 #ifdef MPI
-      nscp_int_tot=(npept-iscp+1)*(npept-iscp+1)
-      call int_bounds(nscp_int_tot,my_scp_inds,my_scp_inde)
+      nscp_int_tot=int(npept-iscp+1,8)*int(npept-iscp+1,8)
+      call int_bounds8(nscp_int_tot,my_scp_inds,my_scp_inde)
       if (lprint) write (iout,*) 'Processor',fg_rank,' CG group',kolor,
      &  ' absolute rank',myrank,
      &  ' nscp_int_tot',nscp_int_tot,' my_scp_inds=',my_scp_inds,
@@ -675,20 +677,20 @@ C Partition the SC-p interaction array
       do i=nnt,nct-1
         if (i.lt.nnt+iscp) then
 cd        write (iout,*) 'i.le.nnt+iscp'
-          call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
+          call int_partition8(ind_scpint,my_scp_inds,my_scp_inde,i,
      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,1),
      &      iscpend(i,1),*14)
         else if (i.gt.nct-iscp) then
 cd        write (iout,*) 'i.gt.nct-iscp'
-          call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
+          call int_partition8(ind_scpint,my_scp_inds,my_scp_inde,i,
      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
      &      iscpend(i,1),*14)
         else
-          call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
+          call int_partition8(ind_scpint,my_scp_inds,my_scp_inde,i,
      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
      &      iscpend(i,1),*14)
           ii=nscp_gr(i)+1
-          call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
+          call int_partition8(ind_scpint,my_scp_inds,my_scp_inde,i,
      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,ii),
      &      iscpend(i,ii),*14)
         endif
@@ -719,7 +721,7 @@ cd        write (iout,*) 'i.gt.nct-iscp'
       if (lprint) then
         write (iout,'(a)') 'SC-p interaction array:'
         do i=iatscp_s,iatscp_e
-          write (iout,'(i3,2(2x,2i3))') 
+          write (iout,'(i7,2(2x,2i7))') 
      &         i,(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i))
         enddo
       endif ! lprint
@@ -782,24 +784,24 @@ c     &  " ivec_start",ivec_start," ivec_end",ivec_end
       endif
 c      nsumgrad=(nres-nnt)*(nres-nnt+1)/2
 c      nlen=nres-nnt+1
-      nsumgrad=(nres-nnt)*(nres-nnt+1)/2
-      nlen=nres-nnt+1
-      call int_bounds(nsumgrad,ngrad_start,ngrad_end)
-      igrad_start=((2*nlen+1)
-     &    -sqrt(float((2*nlen-1)**2-8*(ngrad_start-1))))/2
-      jgrad_start(igrad_start)=
-     &    ngrad_start-(2*nlen-igrad_start)*(igrad_start-1)/2
-     &    +igrad_start
-      jgrad_end(igrad_start)=nres
-      igrad_end=((2*nlen+1)
-     &    -sqrt(float((2*nlen-1)**2-8*(ngrad_end-1))))/2
-      if (igrad_end.gt.igrad_start) jgrad_start(igrad_end)=igrad_end+1
-      jgrad_end(igrad_end)=ngrad_end-(2*nlen-igrad_end)*(igrad_end-1)/2
-     &    +igrad_end
-      do i=igrad_start+1,igrad_end-1
-        jgrad_start(i)=i+1
-        jgrad_end(i)=nres
-      enddo
+c      nsumgrad=(nres-nnt)*(nres-nnt+1)/2
+c      nlen=nres-nnt+1
+c      call int_bounds(nsumgrad,ngrad_start,ngrad_end)
+c      igrad_start=((2*nlen+1)
+c     &    -sqrt(float((2*nlen-1)**2-8*(ngrad_start-1))))/2
+c      jgrad_start(igrad_start)=
+c     &    ngrad_start-(2*nlen-igrad_start)*(igrad_start-1)/2
+c     &    +igrad_start
+c      jgrad_end(igrad_start)=nres
+c      igrad_end=((2*nlen+1)
+c     &    -sqrt(float((2*nlen-1)**2-8*(ngrad_end-1))))/2
+c      if (igrad_end.gt.igrad_start) jgrad_start(igrad_end)=igrad_end+1
+c      jgrad_end(igrad_end)=ngrad_end-(2*nlen-igrad_end)*(igrad_end-1)/2
+c     &    +igrad_end
+c      do i=igrad_start+1,igrad_end-1
+c        jgrad_start(i)=i+1
+c        jgrad_end(i)=nres
+c      enddo
       if (lprint) then 
         write (*,*) 'Processor:',fg_rank,' CG group',kolor,
      & ' absolute rank',myrank,
@@ -818,13 +820,13 @@ c      nlen=nres-nnt+1
      & ' ithetaconstr_start',ithetaconstr_start,' ithetaconstr_end',
      &   ithetaconstr_end
 
-       write (*,*) 'Processor:',fg_rank,myrank,' igrad_start',
-     &   igrad_start,' igrad_end',igrad_end,' ngrad_start',ngrad_start,
-     &   ' ngrad_end',ngrad_end
-       do i=igrad_start,igrad_end
-         write(*,*) 'Processor:',fg_rank,myrank,i,
-     &    jgrad_start(i),jgrad_end(i)
-       enddo
+c       write (*,*) 'Processor:',fg_rank,myrank,' igrad_start',
+c     &   igrad_start,' igrad_end',igrad_end,' ngrad_start',ngrad_start,
+c     &   ' ngrad_end',ngrad_end
+c       do i=igrad_start,igrad_end
+c         write(*,*) 'Processor:',fg_rank,myrank,i,
+c     &    jgrad_start(i),jgrad_end(i)
+c       enddo
       endif
       if (nfgtasks.gt.1) then
         call MPI_Allgather(ivec_start,1,MPI_INTEGER,ivec_displ(0),1,
@@ -1502,6 +1504,31 @@ c---------------------------------------------------------------------------
       return
       end
 c---------------------------------------------------------------------------
+      subroutine int_bounds8(total_ints,lower_bound,upper_bound)
+      implicit none
+      include 'DIMENSIONS'
+      include 'mpif.h'
+      include 'COMMON.SETUP'
+      integer*8 total_ints,lower_bound,upper_bound,nint
+      integer*8 int4proc(0:max_fg_procs),sint4proc(0:max_fg_procs)
+      integer i,nexcess
+      nint=total_ints/nfgtasks
+      do i=1,nfgtasks
+        int4proc(i-1)=nint
+      enddo
+      nexcess=total_ints-nint*nfgtasks
+      do i=1,nexcess
+        int4proc(nfgtasks-i)=int4proc(nfgtasks-i)+1
+      enddo
+      lower_bound=0
+      do i=0,fg_rank-1
+        lower_bound=lower_bound+int4proc(i)
+      enddo 
+      upper_bound=lower_bound+int4proc(fg_rank)
+      lower_bound=lower_bound+1
+      return
+      end
+c---------------------------------------------------------------------------
       subroutine int_bounds1(total_ints,lower_bound,upper_bound)
       implicit none
       include 'DIMENSIONS'
@@ -1566,6 +1593,47 @@ c---------------------------------------------------------------------------
       endif
       return
       end
+c---------------------------------------------------------------------------
+      subroutine int_partition8(int_index,lower_index,upper_index,atom,
+     & at_start,at_end,first_atom,last_atom,int_gr,jat_start,jat_end,*)
+      implicit none
+      include 'DIMENSIONS'
+      include 'COMMON.IOUNITS'
+      integer*8 int_index,lower_index,upper_index
+      integer atom,at_start,at_end,
+     & first_atom,last_atom,int_gr,jat_start,jat_end,int_index_old
+      logical lprn
+      lprn=.false.
+      if (lprn) write (iout,*) 'int_index=',int_index
+      int_index_old=int_index
+      int_index=int_index+last_atom-first_atom+1
+      if (lprn) 
+     &   write (iout,*) 'int_index=',int_index,
+     &               ' int_index_old',int_index_old,
+     &               ' lower_index=',lower_index,
+     &               ' upper_index=',upper_index,
+     &               ' atom=',atom,' first_atom=',first_atom,
+     &               ' last_atom=',last_atom
+      if (int_index.ge.lower_index) then
+        int_gr=int_gr+1
+        if (at_start.eq.0) then
+          at_start=atom
+          jat_start=first_atom-1+lower_index-int_index_old
+        else
+          jat_start=first_atom
+        endif
+        if (lprn) write (iout,*) 'jat_start',jat_start
+        if (int_index.ge.upper_index) then
+          at_end=atom
+          jat_end=first_atom-1+upper_index-int_index_old
+          return1
+        else
+          jat_end=last_atom
+        endif
+        if (lprn) write (iout,*) 'jat_end',jat_end
+      endif
+      return
+      end
 #endif
 c------------------------------------------------------------------------------
       subroutine hpb_partition
index 1180645..6dd113b 100644 (file)
@@ -8,7 +8,7 @@ c-------------------------------------------------------------------------
        include 'DIMENSIONS'
 #ifdef MPI
        include 'mpif.h'
-       integer time00
+       double precision time00
 #endif
        include 'COMMON.VAR'
        include 'COMMON.CHAIN'
@@ -55,7 +55,7 @@ c-------------------------------------------------------------------------
 #endif
 #ifdef FIVEDIAG
       call grad_transform
-      d_a=0.0d0
+      d_a(:,:2*nres)=0.0d0
       if (lprn) then
         write (iout,*) "Potential forces backbone"
         do i=1,nres
@@ -459,9 +459,9 @@ c          write (iout,*) "i",i," itype",itype(i),ntyp1
           endif
         enddo
       enddo
-      DMorig=DM
-      DU1orig=DU1
-      DU2orig=DU2
+      DMorig(:2*nres)=DM(:2*nres)
+      DU1orig(:2*nres)=DU1(:2*nres)
+      DU2orig(:2*nres)=DU2(:2*nres)
       if (gmatout) then
       write (iout,*)"The upper part of the five-diagonal inertia matrix"
       endif
@@ -886,6 +886,9 @@ c---------------------------------------------------------------------------
 c---------------------------------------------------------------------------
       subroutine fivediaginv_mult(ndim,forces,d_a_vec)
       implicit none
+#ifdef MPI
+      include 'mpif.h'
+#endif
       include 'DIMENSIONS'
       include 'COMMON.CHAIN'
       include 'COMMON.IOUNITS'
@@ -896,7 +899,12 @@ c---------------------------------------------------------------------------
       double precision forces(3*ndim),accel(3,0:maxres2),rs(ndim),
      &  xsolv(ndim),d_a_vec(6*nres)
       integer i,j,ind,ichain,n,iposc,innt,inct,inct_prev
-      accel=0.0d0
+#ifdef TIMING
+      include 'COMMON.TIME1'
+      double precision time01
+      time01=MPI_Wtime()
+#endif
+      accel(:,:2*nres)=0.0d0
       do j=1,3
 Compute accelerations in Calpha and SC
         do ichain=1,nchain
@@ -996,6 +1004,9 @@ C Convert d_a to virtual-bon-vector basis
           ind=ind+3
         endif
       enddo
+#ifdef TIMING
+        time_ginvmult=time_ginvmult+MPI_Wtime()-time01
+#endif
 #ifdef DEBUG
       write (iout,*) "d_a_vec"
       write (iout,'(3f10.5)') (d_a_vec(j),j=1,3*(nct-nnt+nside))
index 480aeb2..a1f6b45 100644 (file)
@@ -22,6 +22,7 @@
 !            print *,"START make_SC"
 #ifdef DEBUG
       write (iout,*) "make_SCSC_inter_list maxint_res",maxint_res
+      write (iout,*) "iatsc_s",iatsc_s," iatsc_e",iatsc_e
 #endif
           r_buff_list=5.0d0
             ilist_sc=0
 #ifdef MPI
 #ifdef DEBUG
       write (iout,*) "before MPIREDUCE",ilist_sc
-      do i=1,ilist_sc
-      write (iout,*) i,contlisti(i),contlistj(i)
-      enddo
+c      do i=1,ilist_sc
+c      write (iout,*) i,contlisti(i),contlistj(i)
+c      enddo
 #endif
       if (nfgtasks.gt.1)then
 
@@ -171,10 +172,226 @@ c        write (iout,*) "SCSC after bcast ierr",ierr
       write (iout,*) i,newcontlisti(i),newcontlistj(i)
       enddo
 #endif
-        call int_bounds(g_ilist_sc,g_listscsc_start,g_listscsc_end)
+      call int_bounds(g_ilist_sc,g_listscsc_start,g_listscsc_end)
+#ifdef DEBUG
+      write (iout,*) "g_listscsc_start",g_listscsc_start,
+     &  "g_listscsc_end",g_listscsc_end
       return
+#endif
       end
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+      subroutine make_SCSC_inter_list_RESPA
+      implicit none
+      include "DIMENSIONS"
+#ifdef MPI
+      include 'mpif.h'
+      include "COMMON.SETUP"
+#endif
+      include "COMMON.CONTROL"
+      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_long(maxint_res*maxres),
+     &  contlisti_short(maxint_res*maxres),
+     &  contlistj_long(maxint_res*maxres),
+     &  contlistj_short(maxint_res*maxres)
+!      integer :: newcontlisti(200*nres),newcontlistj(200*nres) 
+      integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint,
+     &  ilist_sc_long,g_ilist_sc_long,ilist_sc_short,g_ilist_sc_short
+      integer displ(0:max_fg_procs),i_ilist_sc_long(0:max_fg_procs),
+     & i_ilist_sc_short(0:max_fg_procs),ierr
+      logical lprn /.false./
+      double precision boxshift
+      double precision d_scale,r_respa_buf
+!            print *,"START make_SC"
+#ifdef DEBUG
+      write (iout,*) "make_SCSC_inter_list maxint_res",maxint_res
+      write (iout,*) "iatsc_s",iatsc_s," iatsc_e",iatsc_e
+#endif
+      r_buff_list=5.0d0
+      r_respa_buf=rlamb
+      ilist_sc_long=0
+      ilist_sc_short=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)
+        call to_box(xi,yi,zi)
+        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)
+            call to_box(xj,yj,zj)
+            xj=boxshift(xj-xi,boxxsize)
+            yj=boxshift(yj-yi,boxysize)
+            zj=boxshift(zj-zi,boxzsize)
+            dist_init=dsqrt(xj*xj+yj*yj+zj*zj)
+! r_buff_list is a read value for a buffer 
+            if (dist_init.le.(r_cut_int+r_buff_list)) then
+! Here the list is created
+              d_scale=dist_init/sigmaii(itypi,itypj)
+              if (d_scale.le.r_cut_respa+r_respa_buf) then
+                ilist_sc_short=ilist_sc_short+1
+                contlisti_short(ilist_sc_short)=i
+                contlistj_short(ilist_sc_short)=j
+              endif
+              if (d_scale.gt.r_cut_respa-rlamb-r_respa_buf) then
+                ilist_sc_long=ilist_sc_long+1
+! this can be substituted by cantor and anti-cantor
+                contlisti_long(ilist_sc_long)=i
+                contlistj_long(ilist_sc_long)=j
+              endif
+            endif
+          enddo
+        enddo
+      enddo
+#ifdef MPI
+#ifdef DEBUG
+      write (iout,*) "before MPIREDUCE ilist_sc_long",ilist_sc_long
+c      do i=1,ilist_sc_long
+c      write (iout,*) i,contlisti_long(i),contlistj_long(i)
+c      enddo
+      write (iout,*) "before MPIREDUCE ilist_sc_short",ilist_sc_short
+c      do i=1,ilist_sc_short
+c      write (iout,*) i,contlisti_short(i),contlistj_short(i)
+c      enddo
+#endif
+      if (nfgtasks.gt.1)then
+
+        call MPI_Reduce(ilist_sc_long,g_ilist_sc_long,1,
+     &                  MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
+        call MPI_Reduce(ilist_sc_short,g_ilist_sc_short,1,
+     &                  MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
+c        write (iout,*) "SCSC after reduce ierr",ierr
+        if (fg_rank.eq.0.and.(g_ilist_sc_long.gt.maxres*maxint_res .or. 
+     &      g_ilist_sc_short.gt.maxres*maxint_res)) then
+          if ((me.eq.king.or.out1file).and.energy_dec) then
+            write (iout,*) "Too many SCSC interactions",
+     &      g_ilist_sc_long,g_ilist_sc_short,
+     &       " only",maxres*maxint_res," allowed."
+            write (iout,*) "Specify a smaller r_cut_int and resubmit"
+            call flush(iout)
+          endif
+          write (*,*) "Processor:",me,": Too many SCSC interactions",
+     &      g_ilist_sc_long+g_ilist_sc_short," only",
+     &      maxres*maxint_res," allowed."
+            write (*,*) "Specify a smaller r_cut_int and resubmit"
+          call MPI_Abort(MPI_COMM_WORLD,ierr)
+        endif
+c        write(iout,*) "before bcast",g_ilist_sc_long
+        call MPI_Gather(ilist_sc_long,1,MPI_INTEGER,
+     &                  i_ilist_sc_long,1,MPI_INTEGER,king,FG_COMM,IERR)
+c        write (iout,*) "SCSC after gather ierr",ierr
+        displ(0)=0
+        do i=1,nfgtasks-1,1
+          displ(i)=i_ilist_sc_long(i-1)+displ(i-1)
+        enddo
+!        write(iout,*) "before gather",displ(0),displ(1)        
+        call MPI_Gatherv(contlisti_long,ilist_sc_long,MPI_INTEGER,
+     &             newcontlisti_long,i_ilist_sc_long,displ,MPI_INTEGER,
+     &             king,FG_COMM,IERR)
+c        write (iout,*) "SCSC after gatherv ierr",ierr
+        call MPI_Gatherv(contlistj_long,ilist_sc_long,MPI_INTEGER,
+     &             newcontlistj_long,i_ilist_sc_long,displ,MPI_INTEGER,
+     &             king,FG_COMM,IERR)
+        call MPI_Bcast(g_ilist_sc_long,1,MPI_INT,king,FG_COMM,IERR)
+c        write (iout,*) "SCSC bcast reduce ierr",ierr
+!        write(iout,*) "before bcast",g_ilist_sc_long
+!        call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+        call MPI_Bcast(newcontlisti_long,g_ilist_sc_long,MPI_INT,king,
+     &       FG_COMM,IERR)
+c        write (iout,*) "SCSC bcast reduce ierr",ierr
+        call MPI_Bcast(newcontlistj_long,g_ilist_sc_long,MPI_INT,king,
+     &       FG_COMM,IERR)
+c        write (iout,*) "SCSC after bcast ierr",ierr
+!        write(iout,*) "before gather",displ(0),displ(1)        
+c        write(iout,*) "before bcast",g_ilist_sc_short
+        call MPI_Gather(ilist_sc_short,1,MPI_INTEGER,
+     &                i_ilist_sc_short,1,MPI_INTEGER,king,FG_COMM,IERR)
+c        write (iout,*) "SCSC after gather ierr",ierr
+        displ(0)=0
+        do i=1,nfgtasks-1,1
+          displ(i)=i_ilist_sc_short(i-1)+displ(i-1)
+        enddo
+!        write(iout,*) "before gather",displ(0),displ(1)        
+        call MPI_Gatherv(contlisti_short,ilist_sc_short,MPI_INTEGER,
+     &            newcontlisti_short,i_ilist_sc_short,displ,MPI_INTEGER,
+     &            king,FG_COMM,IERR)
+c        write (iout,*) "SCSC after gatherv ierr",ierr
+        call MPI_Gatherv(contlistj_short,ilist_sc_short,MPI_INTEGER,
+     &           newcontlistj_short,i_ilist_sc_short,displ,MPI_INTEGER,
+     &           king,FG_COMM,IERR)
+        call MPI_Bcast(g_ilist_sc_short,1,MPI_INT,king,FG_COMM,IERR)
+c        write (iout,*) "SCSC bcast reduce ierr",ierr
+!        write(iout,*) "before bcast",g_ilist_sc_short
+!        call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+        call MPI_Bcast(newcontlisti_short,g_ilist_sc_short,MPI_INT,king,
+     &       FG_COMM,IERR)
+c        write (iout,*) "SCSC bcast reduce ierr",ierr
+        call MPI_Bcast(newcontlistj_short,g_ilist_sc_short,MPI_INT,king,
+     ^       FG_COMM,IERR)
+c        write (iout,*) "SCSC after bcast ierr",ierr
+        else
+#endif
+          g_ilist_sc_long=ilist_sc_long
+
+          do i=1,ilist_sc_long
+            newcontlisti_long(i)=contlisti_long(i)
+            newcontlistj_long(i)=contlistj_long(i)
+          enddo
+
+          g_ilist_sc_short=ilist_sc_short
+
+          do i=1,ilist_sc_short
+            newcontlisti_short(i)=contlisti_short(i)
+            newcontlistj_short(i)=contlistj_short(i)
+          enddo
+#ifdef MPI
+        endif
+#endif      
+      if (fg_rank.eq.0 .and. (me.eq.king.or.out1file) .and. energy_dec) 
+     & write (iout,'(a30,2i10,a,2i4)') 
+     &  "Number of long- and short-range SC-SC interactions",
+     &  g_ilist_sc_long,g_ilist_sc_short," per residue on average",
+     &  g_ilist_sc_long/nres,g_ilist_sc_short/nres
+#ifdef DEBUG
+      write (iout,*) 
+     &  "make_SCSC_inter_list: g_ilist_sc_long after GATHERV",
+     &  g_ilist_sc_long
+      write (iout,*) "List of long-range SCSC interactions"
+      do i=1,g_ilist_sc_long
+      write (iout,*) i,newcontlisti_long(i),newcontlistj_long(i)
+      enddo
+      write (iout,*) 
+     &  "make_SCSC_inter_list: g_ilist_sc_short after GATHERV",
+     &  g_ilist_sc_short
+      write (iout,*) "List of short-range SCSC interactions"
+      do i=1,g_ilist_sc_short
+      write (iout,*) i,newcontlisti_short(i),newcontlistj_short(i)
+      enddo
+#endif
+      call int_bounds(g_ilist_sc_long,g_listscsc_start_long,
+     & g_listscsc_end_long)
+      call int_bounds(g_ilist_sc_short,g_listscsc_start_short,
+     & g_listscsc_end_short)
+#ifdef DEBUG
+      write (iout,*) "g_list_sc_start",g_listscsc_start_long,
+     &  "g_list_sc_end",g_listscsc_end_long
+      write (iout,*)"g_list_sc_start_short",g_listscsc_start_short,
+     &  "g_list_sc_end_short",g_listscsc_end_short
+#endif
+      return
+      end
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       subroutine make_SCp_inter_list
       implicit none
       include "DIMENSIONS"
@@ -204,8 +421,8 @@ c     &  contlistscpj_f(2*maxint_res*maxres)
       write (iout,*) "make_SCp_inter_list maxint_res",maxint_res
 #endif
       r_buff_list=5.0
-            ilist_scp=0
-            ilist_scp_first=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))
@@ -377,11 +594,234 @@ c        write (iout,*) "SCp bcast reduce ierr",ierr
 !      enddo
 #endif
         call int_bounds(g_ilist_scp,g_listscp_start,g_listscp_end)
+#ifdef DEBUG
+      write (iout,*) "g_listscp_start",g_listscp_start,
+     &  "g_listscp_end",g_listscp_end
+#endif
+      return
+      end 
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+      subroutine make_SCp_inter_list_RESPA
+      implicit none
+      include "DIMENSIONS"
+#ifdef MPI
+      include 'mpif.h'
+      include "COMMON.SETUP"
+#endif
+      include "COMMON.CONTROL"
+      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_long(2*maxint_res*maxres),
+     & contlistscpi_short(2*maxint_res*maxres),
+     & contlistscpj_long(2*maxint_res*maxres),
+     & contlistscpj_short(2*maxint_res*maxres)
+!      integer :: newcontlistscpi(200*nres),newcontlistscpj(200*nres)
+      integer i,j,iteli,itypj,subchap,xshift,yshift,zshift,iint,
+     & ilist_scp_long,ilist_scp_short,g_ilist_scp_long,g_ilist_scp_short
+      integer displ(0:max_fg_procs),i_ilist_scp_long(0:max_fg_procs),
+     & i_ilist_scp_short(0:max_fg_procs),ierr
+c      integer contlistscpi_f(2*maxint_res*maxres),
+c     &  contlistscpj_f(2*maxint_res*maxres)
+      double precision boxshift
+      double precision d_scale,r_respa_buf
+!            print *,"START make_SC"
+#ifdef DEBUG
+      write (iout,*) "make_SCp_inter_list maxint_res",maxint_res
+#endif
+      r_buff_list=5.0
+      r_respa_buf=rlamb
+      ilist_scp_long=0
+      ilist_scp_short=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))
+        call to_box(xi,yi,zi)
+        iteli=itel(i)
+        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)
+            call to_box(xj,yj,zj)
+            xj=boxshift(xj-xi,boxxsize)
+            yj=boxshift(yj-yi,boxysize)
+            zj=boxshift(zj-zi,boxzsize)
+            dist_init=dsqrt(xj*xj+yj*yj+zj*zj)
+! r_buff_list is a read value for a buffer 
+            if (dist_init.le.(r_cut_int+r_buff_list)) then
 
+              d_scale=dist_init/rscp(itypj,iteli)
+              if (d_scale.le.r_cut_respa+r_respa_buf) then
+! Here the list is created
+                ilist_scp_short=ilist_scp_short+1
+                contlistscpi_short(ilist_scp_short)=i
+                contlistscpj_short(ilist_scp_short)=j
+              endif
+              if (d_scale.gt.r_cut_respa-rlamb-r_respa_buf) then
+! this can be substituted by cantor and anti-cantor
+                ilist_scp_long=ilist_scp_long+1
+                contlistscpi_long(ilist_scp_long)=i
+                contlistscpj_long(ilist_scp_long)=j
+              endif
+            endif
+          enddo
+        enddo
+      enddo
+#ifdef MPI
+#ifdef DEBUG
+      write (iout,*) "before MPIREDUCE",ilist_scp_long,ilist_scp_short
+      write (iout,*) "Long-range scp interaction list"
+      do i=1,ilist_scp_long
+        write (iout,*) i,contlistscpi_long(i),contlistscpj_long(i)
+      enddo
+      write (iout,*) "Short-range scp interaction list"
+      do i=1,ilist_scp_short
+        write (iout,*) i,contlistscpi_short(i),contlistscpj_short(i)
+      enddo
+#endif
+      if (nfgtasks.gt.1)then
+
+        call MPI_Reduce(ilist_scp_long,g_ilist_scp_long,1,
+     &    MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
+        call MPI_Reduce(ilist_scp_short,g_ilist_scp_short,1,
+     &    MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
+c        write (iout,*) "SCp after reduce ierr",ierr
+        if (fg_rank.eq.0.and.(g_ilist_scp_long.gt.
+     &      2*maxres*maxint_res .or. g_ilist_scp_short.gt.
+     &      2*maxres*maxint_res)) then
+          if ((me.eq.king.or.out1file).and.energy_dec) then
+            write (iout,*) "Too many SCp interactions",
+     &      g_ilist_scp_long+g_ilist_scp_short," only",
+     &      2*maxres*maxint_res," allowed."
+            write (iout,*) "Specify a smaller r_cut_int and resubmit"
+            call flush(iout)
+          endif
+          write (*,*) "Processor:",me,": Too many SCp interactions",
+     &      g_ilist_scp_long+g_ilist_scp_short," only",
+     &      2*maxres*maxint_res," allowed."
+          write (*,*) "Specify a smaller r_cut_int and resubmit"
+          call MPI_Abort(MPI_COMM_WORLD,ierr)
+        endif
+c        write(iout,*) "before bcast",g_ilist_sc
+        call MPI_Gather(ilist_scp_long,1,MPI_INTEGER,
+     &               i_ilist_scp_long,1,MPI_INTEGER,king,FG_COMM,IERR)
+        call MPI_Gather(ilist_scp_short,1,MPI_INTEGER,
+     &               i_ilist_scp_short,1,MPI_INTEGER,king,FG_COMM,IERR)
+c        write (iout,*) "SCp after gather ierr",ierr
+        displ(0)=0
+        do i=1,nfgtasks-1,1
+          displ(i)=i_ilist_scp_long(i-1)+displ(i-1)
+        enddo
+!        write(iout,*) "before gather",displ(0),displ(1)
+        call MPI_Gatherv(contlistscpi_long,ilist_scp_long,MPI_INTEGER,
+     &         newcontlistscpi_long,i_ilist_scp_long,displ,MPI_INTEGER,
+     &         king,FG_COMM,IERR)
+c        write (iout,*) "SCp after gatherv ierr",ierr
+        call MPI_Gatherv(contlistscpj_long,ilist_scp_long,MPI_INTEGER,
+     &         newcontlistscpj_long,i_ilist_scp_long,displ,MPI_INTEGER,
+     &         king,FG_COMM,IERR)
+c        write (iout,*) "SCp after gatherv ierr",ierr
+        call MPI_Bcast(g_ilist_scp_long,1,MPI_INT,king,FG_COMM,IERR)
+c        write (iout,*) "SCp after bcast ierr",ierr
+!        write(iout,*) "before bcast",g_ilist_sc
+!        call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+        call MPI_Bcast(newcontlistscpi_long,g_ilist_scp_long,MPI_INT,
+     &                   king,FG_COMM,IERR)
+c        write (iout,*) "SCp after bcast ierr",ierr
+        call MPI_Bcast(newcontlistscpj_long,g_ilist_scp_long,MPI_INT,
+     &                   king,FG_COMM,IERR)
+c        write (iout,*) "SCp bcast reduce ierr",ierr
+!        call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+        displ(0)=0
+        do i=1,nfgtasks-1,1
+          displ(i)=i_ilist_scp_short(i-1)+displ(i-1)
+        enddo
+!        write(iout,*) "before gather",displ(0),displ(1)
+        call MPI_Gatherv(contlistscpi_short,ilist_scp_short,MPI_INTEGER,
+     &        newcontlistscpi_short,i_ilist_scp_short,displ,MPI_INTEGER,
+     &        king,FG_COMM,IERR)
+c        write (iout,*) "SCp after gatherv ierr",ierr
+        call MPI_Gatherv(contlistscpj_short,ilist_scp_short,MPI_INTEGER,
+     &        newcontlistscpj_short,i_ilist_scp_short,displ,MPI_INTEGER,
+     &        king,FG_COMM,IERR)
+c        write (iout,*) "SCp after gatherv ierr",ierr
+        call MPI_Bcast(g_ilist_scp_short,1,MPI_INT,king,FG_COMM,IERR)
+        call MPI_Bcast(newcontlistscpi_short,g_ilist_scp_short,MPI_INT,
+     &        king,FG_COMM,IERR)
+c        write (iout,*) "SCp after bcast ierr",ierr
+        call MPI_Bcast(newcontlistscpj_short,g_ilist_scp_short,MPI_INT,
+     &        king,FG_COMM,IERR)
+      else
+#endif
+        g_ilist_scp_long=ilist_scp_long
+
+        do i=1,ilist_scp_long
+          newcontlistscpi_long(i)=contlistscpi_long(i)
+          newcontlistscpj_long(i)=contlistscpj_long(i)
+        enddo
+        g_ilist_scp_short=ilist_scp_short
+
+        do i=1,ilist_scp_short
+          newcontlistscpi_short(i)=contlistscpi_short(i)
+          newcontlistscpj_short(i)=contlistscpj_short(i)
+        enddo
+#ifdef MPI
+      endif
+#endif
+      if (fg_rank.eq.0 .and. (me.eq.king.or.out1file) .and. energy_dec) 
+     &then
+        write (iout,'(a30,i10,a,i4)') 
+     &  "Number of long-range SC-p interactions",
+     &  g_ilist_scp_long," per residue on average",g_ilist_scp_long/nres
+        write (iout,'(a30,i10,a,i4)') 
+     &  "Number of short-range SC-p interactions",
+     &g_ilist_scp_short," per residue on average",g_ilist_scp_short/nres
+      endif
+#ifdef DEBUG
+      write (iout,*) "make_SCp_inter_list: after GATHERV long-range",
+     &   g_ilist_scp_long
+      do i=1,g_ilist_scp_long
+        write (iout,*) i,newcontlistscpi_long(i),newcontlistscpj_long(i)
+      enddo
+      write (iout,*) "make_SCp_inter_list: after GATHERV short-range",
+     &   g_ilist_scp_short
+      do i=1,g_ilist_scp_short
+        write (iout,*) i,newcontlistscpi_short(i),
+     &   newcontlistscpj_short(i)
+      enddo
+#endif
+      call int_bounds(g_ilist_scp_long,g_listscp_start_long,
+     &  g_listscp_end_long)
+      call int_bounds(g_ilist_scp_short,g_listscp_start_short,
+     &  g_listscp_end_short)
+      if (fg_rank.eq.0 .and. (me.eq.king.or.out1file) .and. energy_dec)
+     &then
+        write (iout,*) "g_listscp_start",g_listscp_start_long,
+     &  "g_listscp_end",g_listscp_end_long
+        write (iout,*)"g_listscp_start_short",g_listscp_start_short,
+     &  "g_listscp_end_short",g_listscp_end_short
+      endif
       return
       end 
 !-----------------------------------------------------------------------------
-      subroutine make_pp_vdw_inter_list
+      subroutine make_pp_vdw_inter_list_RESPA
       implicit none
       include "DIMENSIONS"
 #ifdef MPI
@@ -398,164 +838,148 @@ c        write (iout,*) "SCp bcast reduce ierr",ierr
       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 contlistpp_vdwi(maxint_res*maxres),
-     & contlistpp_vdwj(maxint_res*maxres)
-!      integer :: newcontlistppi(200*nres),newcontlistppj(200*nres)
+      double precision dxj,dyj,dzj
+      integer contlistpp_vdwi_short(maxint_res*maxres),
+     & contlistpp_vdwj_short(maxint_res*maxres)
       integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint,
-     &  ilist_pp_vdw,g_ilist_pp_vdw
-      integer displ(0:max_fg_procs),i_ilist_pp_vdw(0:max_fg_procs),ierr
+     &  ilist_pp_vdw_short,g_ilist_pp_vdw_short
+      integer displ(0:max_fg_procs),
+     &  i_ilist_pp_vdw_short(0:max_fg_procs),ierr
 !            print *,"START make_SC"
+      double precision boxshift
+      double precision d_scale,r_respa_buf
 #ifdef DEBUG
       write (iout,*) "make_pp_vdw_inter_list"
 #endif
-      ilist_pp_vdw=0
+      ilist_pp_vdw_short=0
       r_buff_list=5.0
+      r_respa_buf=rlamb
       do i=iatel_s_vdw,iatel_e_vdw
         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
+        call to_box(xmedi,ymedi,zmedi)
         do j=ielstart_vdw(i),ielend_vdw(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
+          call to_box(xj,yj,zj)
+          xj=boxshift(xj-xmedi,boxxsize)
+          yj=boxshift(yj-ymedi,boxysize)
+          zj=boxshift(zj-zmedi,boxzsize)
+          dist_init=dsqrt(xj*xj+yj*yj+zj*zj)
 
-          if (dsqrt(dist_init).le.(r_cut_int+r_buff_list)) then
+          if (dist_init.le.(r_cut_int+r_buff_list)) then
+            d_scale=dist_init/rpp(itel(i),itel(j))
+            if (d_scale.le.r_cut_respa+r_respa_buf) then
 ! Here the list is created
-            ilist_pp_vdw=ilist_pp_vdw+1
+              ilist_pp_vdw_short=ilist_pp_vdw_short+1
 ! this can be substituted by cantor and anti-cantor
-            contlistpp_vdwi(ilist_pp_vdw)=i
-            contlistpp_vdwj(ilist_pp_vdw)=j
+              contlistpp_vdwi_short(ilist_pp_vdw_short)=i
+              contlistpp_vdwj_short(ilist_pp_vdw_short)=j
+            endif
           endif
-          enddo
-          enddo
+        enddo
+      enddo
 !             enddo
 #ifdef MPI
 #ifdef DEBUG
-      write (iout,*) "before MPIREDUCE",ilist_pp_vdw
-      do i=1,ilist_pp_vdw
-        write (iout,*) i,contlistpp_vdwi(i),contlistpp_vdwj(i)
+      write (iout,*) "before MPIREDUCE longrange",ilist_pp_vdw_long
+      do i=1,ilist_pp_vdw_long
+        write (iout,*) i,contlistpp_vdwi_long(i),contlistpp_vdwj_long(i)
+      enddo
+      write (iout,*) "before MPIREDUCE shortrange",ilist_pp_vdw_short
+      do i=1,ilist_pp_vdw_short
+        write (iout,*) i,contlistpp_vdwi_short(i),
+     &    contlistpp_vdwj_short(i)
       enddo
 #endif
       if (nfgtasks.gt.1)then
 
-        call MPI_Reduce(ilist_pp_vdw,g_ilist_pp_vdw,1,
+        call MPI_Reduce(ilist_pp_vdw_short,g_ilist_pp_vdw_short,1,
      &    MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
-        if (fg_rank.eq.0.and.g_ilist_pp_vdw.gt.maxres*maxint_res) then
+        if (fg_rank.eq.0.and.g_ilist_pp_vdw_short.gt.maxres*maxint_res) 
+     &  then
           if ((me.eq.king.or.out1file).and.energy_dec) then
             write (iout,*) "Too many pp VDW interactions",
-     &      g_ilist_pp_vdw," only",maxres*maxint_res," allowed."
+     &      g_ilist_pp_vdw_short," only",maxres*maxint_res," allowed."
             write (iout,*) "Specify a smaller r_cut_int and resubmit"
             call flush(iout)
           endif
           write (*,*) "Processor:",me,": Too many pp VDW interactions",
-     &      g_ilist_pp_vdw," only",maxres*maxint_res," allowed."
+     &      g_ilist_pp_vdw_short," only",maxres*maxint_res," allowed."
           write (8,*) "Specify a smaller r_cut_int and resubmit"
           call MPI_Abort(MPI_COMM_WORLD,ierr)
         endif
 !        write(iout,*) "before bcast",g_ilist_sc
-        call MPI_Gather(ilist_pp_vdw,1,MPI_INTEGER,
-     &                  i_ilist_pp_vdw,1,MPI_INTEGER,king,FG_COMM,IERR)
+        call MPI_Gather(ilist_pp_vdw_short,1,MPI_INTEGER,
+     &            i_ilist_pp_vdw_short,1,MPI_INTEGER,king,FG_COMM,IERR)
         displ(0)=0
         do i=1,nfgtasks-1,1
-          displ(i)=i_ilist_pp_vdw(i-1)+displ(i-1)
+          displ(i)=i_ilist_pp_vdw_short(i-1)+displ(i-1)
         enddo
 !        write(iout,*) "before gather",displ(0),displ(1)
-        call MPI_Gatherv(contlistpp_vdwi,ilist_pp_vdw,MPI_INTEGER,
-     &              newcontlistpp_vdwi,i_ilist_pp_vdw,displ,MPI_INTEGER,
-     &              king,FG_COMM,IERR)
-        call MPI_Gatherv(contlistpp_vdwj,ilist_pp_vdw,MPI_INTEGER,
-     &              newcontlistpp_vdwj,i_ilist_pp_vdw,displ,MPI_INTEGER,
-     &              king,FG_COMM,IERR)
-        call MPI_Bcast(g_ilist_pp_vdw,1,MPI_INT,king,FG_COMM,IERR)
+        call MPI_Gatherv(contlistpp_vdwi_short,ilist_pp_vdw_short,
+     &  MPI_INTEGER,newcontlistpp_vdwi_short,i_ilist_pp_vdw_short,displ,
+     &  MPI_INTEGER,king,FG_COMM,IERR)
+        call MPI_Gatherv(contlistpp_vdwj_short,ilist_pp_vdw_short,
+     &  MPI_INTEGER,newcontlistpp_vdwj_short,i_ilist_pp_vdw_short,displ,
+     &  MPI_INTEGER,king,FG_COMM,IERR)
+        call MPI_Bcast(g_ilist_pp_vdw_short,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(newcontlistpp_vdwi,g_ilist_pp_vdw,MPI_INT,king,
-     &                   FG_COMM,IERR)
-        call MPI_Bcast(newcontlistpp_vdwj,g_ilist_pp_vdw,MPI_INT,king,
-     &                   FG_COMM,IERR)
-
+        call MPI_Bcast(newcontlistpp_vdwi_short,g_ilist_pp_vdw_short,
+     &   MPI_INT,king,FG_COMM,IERR)
+        call MPI_Bcast(newcontlistpp_vdwj_short,g_ilist_pp_vdw_short,
+     &   MPI_INT,king,FG_COMM,IERR)
 !        call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
-
-        else
+      else
 #endif
-        g_ilist_pp_vdw=ilist_pp_vdw
+        g_ilist_pp_vdw_short=ilist_pp_vdw_short
 
-        do i=1,ilist_pp_vdw
-          newcontlistpp_vdwi(i)=contlistpp_vdwi(i)
-          newcontlistpp_vdwj(i)=contlistpp_vdwj(i)
+        do i=1,ilist_pp_vdw_short
+          newcontlistpp_vdwi_short(i)=contlistpp_vdwi_short(i)
+          newcontlistpp_vdwj_short(i)=contlistpp_vdwj_short(i)
         enddo
 #ifdef MPI
-        endif
+      endif
 #endif
-        call int_bounds(g_ilist_pp_vdw,g_listpp_vdw_start,
-     &       g_listpp_vdw_end)
-      if (fg_rank.eq.0 .and. (me.eq.king.or.out1file) .and. energy_dec) 
-     &write (iout,'(a30,i10,a,i4)') "Number of p-p VDW interactions",
-     & g_ilist_pp_vdw," per residue on average",g_ilist_pp_vdw/nres
+      if (fg_rank.eq.0 .and. (me.eq.king.or.out1file) .and. energy_dec)
+     &then
+      write (iout,*) "Number of short-range p-p VDW interactions",
+     & g_ilist_pp_vdw_short," per residue on average",
+     & g_ilist_pp_vdw_short/nres
+      endif
 #ifdef DEBUG
-      write (iout,*) "g_listpp_vdw_start",g_listpp_vdw_start,
-     &  "g_listpp_vdw_end",g_listpp_vdw_end
+      write (iout,*) "Short-range pp_vdw"
       write (iout,*) "make_pp_vdw_inter_list: after GATHERV",
-     &  g_ilist_pp_vdw
-      do i=1,g_ilist_pp_vdw
-        write (iout,*) i,newcontlistpp_vdwi(i),newcontlistpp_vdwj(i)
+     &  g_ilist_pp_vdw_short
+      do i=1,g_ilist_pp_vdw_short
+        write (iout,*) i,newcontlistpp_vdwi_short(i),
+     &     newcontlistpp_vdwj_short(i)
       enddo
 #endif
+      call int_bounds(g_ilist_pp_vdw_short,g_listpp_vdw_start_short,
+     &       g_listpp_vdw_end_short)
+      if (fg_rank.eq.0 .and. (me.eq.king.or.out1file) .and. energy_dec)
+     &then
+        write (iout,*)"g_listpp_vdw_start_short",
+     &  g_listpp_vdw_start_short,
+     &  "g_listpp_vdw_end_short",g_listpp_vdw_end_short
+      endif
       return
       end
 !-----------------------------------------------------------------------------
@@ -716,18 +1140,17 @@ c       write (iout,*) "After bcast ierr",ierr
 
 !        call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
 
-        else
+      else
 #endif
         g_ilist_pp=ilist_pp
 
         do i=1,ilist_pp
-        newcontlistppi(i)=contlistppi(i)
-        newcontlistppj(i)=contlistppj(i)
+          newcontlistppi(i)=contlistppi(i)
+          newcontlistppj(i)=contlistppj(i)
         enddo
 #ifdef MPI
-        endif
+      endif
 #endif
-        call int_bounds(g_ilist_pp,g_listpp_start,g_listpp_end)
       if (fg_rank.eq.0 .and. (me.eq.king.or.out1file) .and. energy_dec) 
      & write (iout,'(a30,i10,a,i4)') "Number of p-p interactions",
      & g_ilist_pp," per residue on average",g_ilist_pp/nres
@@ -737,5 +1160,11 @@ c       write (iout,*) "After bcast ierr",ierr
       write (iout,*) i,newcontlistppi(i),newcontlistppj(i)
       enddo
 #endif
+      call int_bounds(g_ilist_pp,g_listpp_start,g_listpp_end)
+      if (fg_rank.eq.0 .and. (me.eq.king.or.out1file) .and. energy_dec)
+     &then
+        write (iout,*) "g_listpp_start",g_listpp_start,
+     &  "g_listpp_end",g_listpp_end
+      endif
       return
       end
index 41a1a27..0b16d55 100644 (file)
@@ -197,11 +197,11 @@ c----------------------------------------------------------------------------
 #endif
       include 'COMMON.TIME1'
       double precision z(maxres6),d_a_tmp(maxres6)
-      double precision edum(0:n_ene),time_order(0:10)
+      double precision edum(0:n_ene),time_order(0:11)
 c      double precision Gcopy(maxres2,maxres2)
 c      common /przechowalnia/ Gcopy
       integer icall /0/
-      integer i,j,iorder
+      integer i,j,iorder,ioverlap(maxres),ioverlap_last
 C Workers wait for variables and NF, and NFL from the boss 
       iorder=0
       do while (iorder.ge.0)
@@ -301,6 +301,8 @@ c           call MATOUT2(my_ng_count,dimen3,maxres2,maxers2,ginv_block)
 #endif
         else if (iorder.eq.10) then
           call setup_fricmat
+        else if (iorder.eq.11) then
+          call overlap_sc_list(ioverlap,ioverlap_last,.false.)
         endif
       enddo
       write (*,*) 'Processor',fg_rank,' CG group',kolor,
diff --git a/source/unres/src-HCD-5D/orig_frame_chain.F b/source/unres/src-HCD-5D/orig_frame_chain.F
new file mode 100644 (file)
index 0000000..07053ce
--- /dev/null
@@ -0,0 +1,85 @@
+      subroutine orig_frame_chain(istart)
+C
+C Define the origin and orientation of the coordinate system starting
+C at residue istart and locate sites istart+1 and istart+2. The
+C coordinates of site istart and the respective dc and dc_norm must be
+C pre-defined
+C
+      implicit none
+      integer i,j,istart,ichain
+      double precision cost,sint,cosg,sing,aux
+      include 'DIMENSIONS'
+      include 'COMMON.CHAIN'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.LOCAL'
+      include 'COMMON.GEO'
+      include 'COMMON.VAR'
+      cost=dc_norm(1,istart)
+      aux=dsqrt(dc_norm(2,istart)**2+dc_norm(3,istart)**2)
+      cosg=dc_norm(2,istart)/aux
+      sing=dc_norm(3,istart)/aux
+      cost=dcos(theta(istart+2))
+      sint=dsin(theta(istart+2))
+      prod(1,1,istart)=cost
+      prod(1,2,istart)=-sint
+      prod(1,3,istart)=0.0d0
+      prod(2,1,istart)=sint*cosg
+      prod(2,2,istart)=cost*cosg
+      prod(2,3,istart)=-sing
+      prod(3,1,istart)=sint*sing
+      prod(3,2,istart)=cost*sing
+      prod(3,3,istart)=cosg
+      aux=prod(1,1,istart)*(prod(2,2,istart)*prod(3,3,istart)
+     & -prod(3,2,istart)*prod(2,3,istart))
+     & -prod(1,2,istart)*(prod(2,1,istart)*prod(3,3,istart)
+     & -prod(3,1,istart)*prod(2,3,istart))
+     & +prod(1,3,istart)*(prod(2,1,istart)*prod(3,2,istart)
+     & -prod(3,1,istart)*prod(2,2,istart))
+c      write (iout,*) "orig_frame_chain prod",istart
+c      do i=1,3
+c        write(iout,'(i5,3f10.5)') i,(prod(i,j,istart),j=1,3)
+c      enddo
+c      write (iout,*) "orig_frame_chain: prod",istart," determinant",aux
+      t(1,1,istart)=cost
+      t(1,2,istart)=-sint 
+      t(1,3,istart)= 0.0D0
+      t(2,1,istart)=sint
+      t(2,2,istart)=cost
+      t(2,3,istart)= 0.0D0
+      t(3,1,istart)= 0.0D0
+      t(3,2,istart)= 0.0D0
+      t(3,3,istart)= 1.0D0
+      r(1,1,istart)= 1.0D0
+      r(1,2,istart)= 0.0D0
+      r(1,3,istart)= 0.0D0
+      r(2,1,istart)= 0.0D0
+      r(2,2,istart)= 1.0D0
+      r(2,3,istart)= 0.0D0
+      r(3,1,istart)= 0.0D0
+      r(3,2,istart)= 0.0D0
+      r(3,3,istart)= 1.0D0
+      do i=1,3
+        do j=1,3
+          rt(i,j,istart)=t(i,j,istart)
+        enddo
+      enddo
+      call matmult(prod(1,1,istart),rt(1,1,istart),prod(1,1,istart+1))
+c      aux=prod(1,1,istart+1)*(prod(2,2,istart+1)*prod(3,3,istart+1)
+c     & -prod(3,2,istart+1)*prod(2,3,istart+1))
+c     & -prod(1,2,istart+1)*(prod(2,1,istart+1)*prod(3,3,istart+1)
+c     & -prod(3,1,istart+1)*prod(2,3,istart+1))
+c     & +prod(1,3,istart+1)*(prod(2,1,istart+1)*prod(3,2,istart+1)
+c     & -prod(3,1,istart+1)*prod(2,2,istart+1))
+c      write (iout,*) "orig_frame_chain prod",istart+1
+c      do i=1,3
+c        write(iout,'(i5,3f10.5)') i,(prod(i,j,istart+1),j=1,3)
+c      enddo
+c      write (iout,*)"orig_frame_chain: prod",istart+1," determinant",aux
+      do j=1,3
+        dc_norm(j,istart+1)=prod(j,1,istart+1)
+        dc(j,istart+1)=vbld(istart+2)*prod(j,1,istart+1)
+        c(j,istart+2)=c(j,istart+1)+dc(j,istart+1)
+      enddo
+      call locate_side_chain(istart+1)
+      return
+      end
index 41fe7f6..12d011e 100644 (file)
@@ -64,9 +64,9 @@ c        call flush(iout)
           iterter(ires_old-1)=1
           itype(ires_old)=ntyp1
           iterter(ires_old)=1
-          ishift1=ishift1+1
+c          ishift1=ishift1+1
           ibeg=2
-          write (iout,*) "Chain ended",ires,ishift,ires_old
+          write (iout,*) "Chain ended",ires,ishift,ires_old,ibeg
           if (unres_pdb) then
             do j=1,3
               dc(j,ires)=sccor(j,iii)
@@ -95,8 +95,8 @@ c          write (iout,*) "! ",atom," !",ires
           read (card(18:20),'(a3)') res
 c          write (iout,*) "ires",ires,ires-ishift+ishift1,
 c     &      " ires_old",ires_old
-c          write (iout,*) "ishift",ishift," ishift1",ishift1
-c          write (iout,*) "IRES",ires-ishift+ishift1,ires_old
+c         write (iout,*) "ishift",ishift," ishift1",ishift1
+c         write (iout,*) "IRES",ires-ishift+ishift1,ires_old
           if (ires-ishift+ishift1.ne.ires_old) then
 ! Calculate the CM of the preceding residue.
 !            if (ibeg.eq.0) call sccenter(ires,iii,sccor)
@@ -115,6 +115,7 @@ c                write (iout,'(i5,3f10.5)') ires,(sccor(j,iii),j=1,3)
               sccalc=.true.
             endif
 ! Start new residue.
+c            write (iout,*) "ibeg",ibeg
             if (res.eq.'Cl-' .or. res.eq.'Na+') then
               ires=ires_old
               cycle
@@ -133,10 +134,13 @@ c              write (iout,*) "BEG ires",ires
             else if (ibeg.eq.2) then
 ! Start a new chain
               ishift=-ires_old+ires-1 !!!!!
-              ishift1=ishift1-1    !!!!!
-c              write (iout,*) "New chain started",ires,ishift,ishift1,"!"
+c              ishift1=ishift1-1    !!!!!
+c              write (iout,*) "New chain started",ires,ires_old,ishift,
+c     &           ishift1
               ires=ires-ishift+ishift1
+              write (iout,*) "New chain started ires",ires
               ires_old=ires
+c              ires=ires_old+1
               ibeg=0
             else
               ishift=ishift-(ires-ishift+ishift1-ires_old-1)
@@ -159,7 +163,9 @@ c          write (2,*) "ires",ires," res ",res!," ity"!,ity
           if (atom.eq.'CA' .or. atom.eq.'CH3' .or. 
      &       res.eq.'NHE'.and.atom(:2).eq.'HN') then
             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
-!            write (iout,*) "backbone ",atom
+            read(card(61:66),*) bfac(ires)
+c            write (iout,*) "backbone ",atom
+c            write (iout,*) ires,res,(c(j,ires),j=1,3)
 #ifdef DEBUG
             write (iout,'(i6,i3,2x,a,3f8.3)') 
      &      ires,itype(ires),res,(c(j,ires),j=1,3)
@@ -181,7 +187,7 @@ c            write (2,*) "iii",iii
         endif
       enddo
    10 if(me.eq.king.or..not.out1file) 
-     & write (iout,'(a,i5)') ' Nres: ',ires
+     & write (iout,'(a,i7)') ' Nres: ',ires
 c      write (iout,*) "iii",iii
 C Calculate dummy residue coordinates inside the "chain" of a multichain
 C system
@@ -229,7 +235,7 @@ C 2/15/2013 by Adam: corrected insertion of the first dummy residue
               e2(3)=0.0d0
             endif
             do j=1,3
-              c(j,i)=c(j,i+1)-1.9d0*e2(j)
+              c(j,i)=c(j,i+1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
             enddo
           else !unres_pdb
            do j=1,3
@@ -316,9 +322,9 @@ C Calculate internal coordinates.
       write (iout,'(/a)')
      &  "Cartesian coordinates of the reference structure"
       write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
-     & "Residue ","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
+     & "Residue   ","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
       do ires=1,nres
-        write (iout,'(a3,1x,i4,3f8.3,5x,3f8.3)')
+        write (iout,'(a3,1x,i6,3f8.3,5x,3f8.3)')
      &    restyp(itype(ires)),ires,(c(j,ires),j=1,3),
      &    (c(j,ires+nres),j=1,3)
       enddo
@@ -412,8 +418,9 @@ C and convert the peptide geometry into virtual-chain geometry.
       character*3 seq,res
       character*5 atom
       character*80 card
-      double precision sccor(3,20)
+      double precision sccor(3,50)
       integer rescode,iterter(maxres)
+      logical zero
       do i=1,maxres
          iterter(i)=0
       enddo
@@ -541,7 +548,7 @@ c            write (iout,*) "sidechain ",atom
         endif
       enddo
    10 if(me.eq.king.or..not.out1file) 
-     & write (iout,'(a,i5)') ' Nres: ',ires
+     & write (iout,'(a,i7)') ' Nres: ',ires
 C Calculate dummy residue coordinates inside the "chain" of a multichain
 C system
       nres=ires
@@ -581,7 +588,7 @@ C 2/15/2013 by Adam: corrected insertion of the first dummy residue
               e2(3)=0.0d0
             endif
             do j=1,3
-              c(j,i)=c(j,i+1)-1.9d0*e2(j)
+              c(j,i)=c(j,i+1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
             enddo
           else !unres_pdb
            do j=1,3
@@ -616,7 +623,7 @@ C 2/15/2013 by Adam: corrected insertion of the last dummy residue
             e2(3)=0.0d0
           endif
           do j=1,3
-            c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
+            c(j,nres)=c(j,nres-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
           enddo
         else
         do j=1,3
@@ -648,7 +655,7 @@ C 2/15/2013 by Adam: corrected insertion of the first dummy residue
             e2(3)=0.0d0
           endif
           do j=1,3
-            c(j,1)=c(j,2)-1.9d0*e2(j)
+            c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/dsqrt(2.0d0)
           enddo
         else
         do j=1,3
@@ -676,6 +683,18 @@ C Calculate internal coordinates.
      &    (c(j,ires+nres),j=1,3)
       enddo
       endif
+      zero=.false.
+      do ires=1,nres
+        zero=zero.or.itype(ires).eq.0
+      enddo
+      if (zero) then
+        write (iout,'(2a)') "Gaps in PDB coordinates detected;",
+     &    " look for ZERO in the control output above."
+        write (iout,'(2a)') "Repair the PDB file using MODELLER",
+     &    " or other softwared and resubmit."
+        call flush(iout)
+        stop
+      endif
 C Calculate internal coordinates.
       call int_from_cart(.true.,out_template_coord)
       call sc_loc_geom(.false.)
@@ -683,6 +702,7 @@ C Calculate internal coordinates.
         thetaref(i)=theta(i)
         phiref(i)=phi(i)
       enddo
+      dc(:,0)=c(:,1)
       do i=1,nres-1
         do j=1,3
           dc(j,i)=c(j,i+1)-c(j,i)
index 4fbc0f1..58b86cc 100644 (file)
@@ -149,6 +149,7 @@ C constrains on theta angles WITH_THETA_CONSTR is the keyword
       out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0
       out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0
       call readi(controlcard,'SYM',symetr,1)
+      call readi(controlcard,'PERMUT',npermut,1)
       call reada(controlcard,'TIMLIM',timlim,2800.0D0) ! default 16 hours
       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
@@ -167,6 +168,8 @@ c       write (iout,'(a,f10.1)') 'Time limit (min):',timlim
 c      endif
       call readi(controlcard,'NZ_START',nz_start,0)
       call readi(controlcard,'NZ_END',nz_end,0)
+      call readi(controlcard,'NRAN_START',nran_start,0)
+      write (iout,*) "nran_start",nran_start
 c      call readi(controlcard,'IZ_SC',iz_sc,0)
       timlim=60.0D0*timlim
       safety = 60.0d0*safety
@@ -481,6 +484,7 @@ C
       call readi(controlcard,"NSTEP",n_timestep,1000000)
       call readi(controlcard,"NTWE",ntwe,100)
       call readi(controlcard,"NTWX",ntwx,1000)
+      call readi(controlcard,"REST_FREQ",irest_freq,1000)
       call reada(controlcard,"DT",d_time,1.0d-1)
       call reada(controlcard,"DVMAX",dvmax,2.0d1)
       call reada(controlcard,"DAMAX",damax,1.0d1)
@@ -556,6 +560,7 @@ c  if performing umbrella sampling, fragments constrained are read from the frag
      &   " A"
        write(iout,'(a60,i5)')"Frequency of updating interaction list",
      &   imatupdate
+       write(iout,'(a60,i5)')"Restart writing frequency",irest_freq
        if (RESPA) then
         write (iout,'(2a,i4,a)') 
      &    "A-MTS algorithm used; initial time step for fast-varying",
@@ -741,7 +746,7 @@ C
       integer ilen
       external ilen
       integer iperm,tperm
-      integer i,j,ii,k,l,itrial,itmp,i1,i2,it1,it2,nres_temp,itemp
+      integer i,j,ii,k,l,itrial,itmp,i1,i2,it1,it2,nres_temp
       double precision sumv
 C
 C Read PDB structure if applicable
@@ -857,16 +862,19 @@ cd      print *,'NNT=',NNT,' NCT=',NCT
      &    chain_border(2,i),chain_border1(1,i),chain_border1(2,i)
       enddo
       call chain_symmetry(nchain,nres,itype,chain_border,
-     &    chain_length,npermchain,tabpermchain)
+     &    chain_length,npermchain,tabpermchain,nchain_group,nequiv,
+     &    iequiv,mapchain)
 c      do i=1,nres
 c        write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain),
 c     &    ii=1,npermchain)
 c      enddo
+#ifdef DEBUG
       write(iout,*) "residue permutations"
       do i=1,nres
         write(iout,*) i,(iperm(i,ii),ii=1,npermchain)
       enddo
       call flush(iout)
+#endif
       if (itype(1).eq.ntyp1) nnt=2
       if (itype(nres).eq.ntyp1) nct=nct-1
       write (iout,*) "nnt",nnt," nct",nct
@@ -971,25 +979,34 @@ C     &   write (iout,*) 'FTORS',ftors
         call reada(weightcard,"WDIHC",wdihc,0.591D0)
         write (iout,*) "Weight of dihedral angle restraints",wdihc
         read(inp,'(9x,3f7.3)') 
+c     &     (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
      &     (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
         write (iout,*) "The secprob array"
         do i=nnt,nct
           write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
         enddo
         ndih_constr=0
+        iconstr_dih=0
         do i=nnt+3,nct
           if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1 
      &    .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
             ndih_constr=ndih_constr+1
             idih_constr(ndih_constr)=i
+            iconstr_dih(i)=ndih_constr
             sumv=0.0d0
             do j=1,3
               vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
               sumv=sumv+vpsipred(j,ndih_constr)
             enddo 
-            do j=1,3
-              vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
-            enddo
+            if (sumv.gt.0.0d0) then
+              do j=1,3
+                vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
+              enddo
+            else
+              vpsipred(1,ndih_constr)=1.0d0
+              vpsipred(2,ndih_constr)=0.0d0
+              vpsipred(3,ndih_constr)=0.0d0
+            endif
             phibound(1,ndih_constr)=phihel*deg2rad
             phibound(2,ndih_constr)=phibet*deg2rad
             sdihed(1,ndih_constr)=sigmahel*deg2rad
@@ -1125,8 +1142,11 @@ c      write (iout,*) "After read_dist_constr nhpb",nhpb
       endif
       write (iout,*) "calling read_saxs_consrtr",nsaxs
       if (nsaxs.gt.0) call read_saxs_constr
-
+c      write (iout,*) "After read_saxs_constr"
+c      call flush(iout) 
       if (constr_homology.gt.0) then
+c        write (iout,*) "Calling read_constr_homology"
+c        call flush(iout)
         call read_constr_homology
         if (indpdb.gt.0 .or. pdbref) then
           do i=1,2*nres
@@ -1184,7 +1204,8 @@ c      write (iout,*) "After read_dist_constr nhpb",nhpb
  335        continue
             unres_pdb=.false.
             nres_temp=nres
-            call readpdb
+c            call readpdb
+            call readpdb_template(nmodel_start+1)
             close(ipdbin)
             if (nres.ge.nres_temp) then
               nmodel_start=nmodel_start+1
@@ -1195,22 +1216,11 @@ c      write (iout,*) "After read_dist_constr nhpb",nhpb
                 enddo
               enddo
             else
-c              itemp=nres
-c              nres=nres_temp
-c              call gen_rand_conf(itemp,*115)
-c              nmodel_start=nmodel_start+1
-c              do i=1,2*nres
-c                do j=1,3
-c                  chomo(j,i,nmodel_start)=c(j,i)
-c                enddo
-c              enddo
-c              goto 116
-  115         if (me.eq.king .or. .not. out1file)
-     &          write (iout,'(a,2i5,1x,a)') 
+              if (me.eq.king .or. .not. out1file)
+     &          write (iout,'(a,2i7,1x,a)') 
      &           "Different number of residues",nres_temp,nres,
      &           " model skipped."
             endif
-  116       continue
             nres=nres_temp
           enddo
   332     continue
@@ -2341,6 +2351,7 @@ c1out       open(iout,file=outname,status='unknown')
 #else
       if (me.eq.king .or. .not.out1file)
      &    open(iout,file=outname,status='unknown')
+#define DEBUG
 #ifdef DEBUG
       if (fg_rank.gt.0) then
         write (liczba,'(i3.3)') myrank/nfgtasks
@@ -2349,6 +2360,7 @@ c1out       open(iout,file=outname,status='unknown')
      &   status='unknown')
       endif
 #endif
+#undef DEBUG
       if(me.eq.king) then
        open(igeom,file=intname,status='unknown',access='append')
        open(ipdb,file=pdbname,status='unknown')
@@ -2734,10 +2746,10 @@ c            write (iout,*) "j",j," k",k
             endif
 #ifdef MPI
             if (.not.out1file .or. me.eq.king) 
-     &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
+     &      write (iout,'(a,3i6,f8.2,1pe12.2)') "+dist.restr ",
      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 #else
-            write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
+            write (iout,'(a,3i6,f8.2,1pe12.2)') "+dist.restr ",
      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 #endif
           enddo
@@ -2781,10 +2793,10 @@ c            write (iout,*) "j",j," k",k
             endif
 #ifdef MPI
             if (.not.out1file .or. me.eq.king)
-     &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
+     &      write (iout,'(a,3i6,f8.2,1pe12.2)') "+dist.restr ",
      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 #else
-            write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
+            write (iout,'(a,3i6,f8.2,1pe12.2)') "+dist.restr ",
      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 #endif
           enddo
@@ -2811,13 +2823,13 @@ c     &    fordepth_peak(nhpb_peak+1),npeak
           ipeak(2,npeak)=i
 #ifdef MPI
           if (.not.out1file .or. me.eq.king)
-     &    write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
+     &    write (iout,'(a,5i6,2f8.2,2f10.5,i5)') "+dist.restr ",
      &     nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
      &     ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
      &     dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
      &     fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
 #else
-          write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
+          write (iout,'(a,5i6,2f8.2,2f10.5,i5)') "+dist.restr ",
      &     nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
      &     ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
      &     dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
@@ -2840,11 +2852,11 @@ c        fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
           irestr_type(nhpb)=11
 #ifdef MPI
           if (.not.out1file .or. me.eq.king)
-     &    write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
+     &    write (iout,'(a,4i6,2f8.2,2f10.5,i5)') "+dist.restr ",
      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
 #else
-          write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
+          write (iout,'(a,4i6,2f8.2,2f10.5,i5)') "+dist.restr ",
      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
 #endif
@@ -2904,12 +2916,12 @@ c     &      dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
           endif
 #ifdef MPI
           if (.not.out1file .or. me.eq.king)
-     &    write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
+     &    write (iout,'(a,4i6,2f8.2,3f10.5,i5)') "+dist.restr ",
      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
      &     irestr_type(nhpb)
 #else
-          write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
+          write (iout,'(a,4i6,2f8.2,3f10.5,i5)') "+dist.restr ",
      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
      &     irestr_type(nhpb)
@@ -2938,10 +2950,10 @@ C        print *,"in else"
         endif
 #ifdef MPI
           if (.not.out1file .or. me.eq.king)
-     &    write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
+     &    write (iout,'(a,4i6,f8.2,f10.1)') "+dist.restr ",
      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
 #else
-          write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
+          write (iout,'(a,4i6,f8.2,f10.1)') "+dist.restr ",
      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
 #endif
         endif
@@ -2984,13 +2996,13 @@ C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
             dhpb(nhpb)=dist(i+nres*ii,j+nres*jj)
 #ifdef MPI
             if (.not.out1file .or. me.eq.king) then
-            write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
+            write (iout,'(a,4i6,2f8.2,3f10.5,i5)') "+dist.restr ",
      &       nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
      &       dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
      &       irestr_type(nhpb)
             endif
 #else
-            write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
+            write (iout,'(a,4i6,2f8.2,3f10.5,i5)') "+dist.restr ",
      &       nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
      &       dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
      &       irestr_type(nhpb)
@@ -3054,8 +3066,8 @@ c    &    sigma_odl_temp(maxres,maxres,max_template)
       character*2 kic2
       character*24 model_ki_dist, model_ki_angle
       character*500 controlcard
-      integer ki,i,ii,j,k,l,ii_in_use(maxdim),i_tmp,idomain_tmp,irec,
-     & ik,iistart,nres_temp
+      integer ki,i,ii,j,k,l,ii_in_use(maxdim_cont),i_tmp,idomain_tmp,
+     & irec,ik,iistart,nres_temp
       integer ilen
       external ilen
       logical liiflag,lfirst
@@ -3108,7 +3120,7 @@ c      endif
          call card_concat(controlcard)
          read(controlcard,*) (waga_homology(i),i=1,homol_nset) 
          if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
-          write(iout,*) "iset homology_weight "
+c          write(iout,*) "iset homology_weight "
           do i=1,homol_nset
            write(iout,*) i,waga_homology(i)
           enddo
@@ -3122,30 +3134,30 @@ c      endif
 cd      write (iout,*) "nnt",nnt," nct",nct
 cd      call flush(iout)
 
+      if (read_homol_frag) then
+        call read_klapaucjusz
+      else
 
       lim_odl=0
       lim_dih=0
 c
 c      write(iout,*) 'nnt=',nnt,'nct=',nct
 c
-      do i = nnt,nct
-        do k=1,constr_homology
-         idomain(k,i)=0
-        enddo
-      enddo
-
-      ii=0
-      do i = nnt,nct-2 
-        do j=i+2,nct 
-        ii=ii+1
-        ii_in_use(ii)=0
-        enddo
-      enddo
-
-      if (read_homol_frag) then
-        call read_klapaucjusz
-      else
+c      do i = nnt,nct
+c        do k=1,constr_homology
+c         idomain(k,i)=0
+c        enddo
+c      enddo
+       idomain=0
 
+c      ii=0
+c      do i = nnt,nct-2 
+c        do j=i+2,nct 
+c        ii=ii+1
+c        ii_in_use(ii)=0
+c        enddo
+c      enddo
+      ii_in_use=0
       do k=1,constr_homology
 
         read(inp,'(a)') pdbfile
@@ -3265,8 +3277,8 @@ c    &                       constr_homology
               endif
               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) 
             else
-              ii=ii+1
-              l_homo(k,ii)=.false.
+c              ii=ii+1
+c              l_homo(k,ii)=.false.
             endif
             enddo
           enddo
@@ -3456,26 +3468,26 @@ cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
        write (iout,*) "Distance restraints from templates"
        do ii=1,lim_odl
-       write(iout,'(3i5,100(2f8.2,1x,l1,4x))') 
+       write(iout,'(3i7,100(2f8.2,1x,l1,4x))') 
      &  ii,ires_homo(ii),jres_homo(ii),
      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
      &  ki=1,constr_homology)
        enddo
        write (iout,*) "Dihedral angle restraints from templates"
        do i=nnt+3,nct
-        write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
+        write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i)),
      &      (rad2deg*dih(ki,i),
      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
        enddo
        write (iout,*) "Virtual-bond angle restraints from templates"
        do i=nnt+2,nct
-        write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
+        write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i)),
      &      (rad2deg*thetatpl(ki,i),
      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
        enddo
        write (iout,*) "SC restraints from templates"
        do i=nnt,nct
-        write(iout,'(i5,100(4f8.2,4x))') i,
+        write(iout,'(i7,100(4f8.2,4x))') i,
      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
        enddo
@@ -3641,14 +3653,17 @@ c----------------------------------------------------------------------
       character*2 kic2
       character*24 model_ki_dist, model_ki_angle
       character*500 controlcard
-      integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp,
-     & ik,ll,ii,kk,iistart,iishift,lim_xx
+      integer ki, i, j, jj,k, l, ii_in_use(maxdim_cont),i_tmp,
+     & idomain_tmp,
+     & ik,ll,lll,ii_old,ii,iii,ichain,kk,iistart,iishift,lim_xx,igr,
+     & i01,i10,nnt_chain,nct_chain
+      integer itype_temp(maxres)
       double precision distal
       logical lprn /.true./
       integer nres_temp
       integer ilen
       external ilen
-      logical liiflag
+      logical liiflag,lfirst
 c
 c
       double precision rescore_tmp,x12,y12,z12,rescore2_tmp
@@ -3662,14 +3677,48 @@ c For new homol impl
 c
       include 'COMMON.VAR'
 c
+c      write (iout,*) "READ_KLAPAUCJUSZ"
+c      print *,"READ_KLAPAUCJUSZ"
+c      call flush(iout)
       call getenv("FRAGFILE",fragfile) 
+      write (iout,*) "Opening", fragfile
+      call flush(iout)
       open(ientin,file=fragfile,status="old",err=10)
-      read(ientin,*) constr_homology,nclust
-      nmodel_start=constr_homology
-      l_homo = .false.
+c      write (iout,*) " opened"
+c      call flush(iout)
+
       sigma_theta=0.0
       sigma_d=0.0
       sigma_dih=0.0
+      l_homo = .false.
+
+      nres_temp=nres
+      itype_temp=itype
+      ii=0
+      lim_odl=0
+
+c      write (iout,*) "Entering loop"
+c      call flush(iout)
+
+      DO IGR = 1,NCHAIN_GROUP
+
+c      write (iout,*) "igr",igr
+      call flush(iout)
+      read(ientin,*) constr_homology,nclust
+
+      if (start_from_model) then
+        nmodel_start=constr_homology
+      else
+        nmodel_start=0
+      endif
+
+      ii_old=lim_odl
+
+      ichain=iequiv(1,igr)
+      nnt_chain=chain_border(1,ichain)-chain_border1(1,ichain)+1
+      nct_chain=chain_border(2,ichain)-chain_border1(1,ichain)+1
+c      write (iout,*) "nnt_chain",nnt_chain," nct_chain",nct_chain
+
 c Read pdb files 
       do k=1,constr_homology 
         read(ientin,'(a)') pdbfile
@@ -3683,10 +3732,10 @@ c Read pdb files
         stop
   34    continue
         unres_pdb=.false.
-        nres_temp=nres
+c        nres_temp=nres
         call readpdb_template(k)
         nres_chomo(k)=nres
-        nres=nres_temp
+c        nres=nres_temp
         do i=1,nres
           rescore(k,i)=0.2d0
           rescore2(k,i)=1.0d0
@@ -3705,6 +3754,7 @@ c
         do ll = 1,ninclust(l)
         
         k = inclust(ll,l)
+c        write (iout,*) "l",l," ll",ll," k",k
         do i=1,nres
           idomain(k,i)=0
         enddo
@@ -3720,7 +3770,7 @@ c     Distance restraints
 c
 c          ... --> odl(k,ii)
 C Copy the coordinates from reference coordinates (?)
-        nres_temp=nres
+c        nres_temp=nres
         nres=nres_chomo(k)
         do i=1,2*nres
           do j=1,3
@@ -3734,11 +3784,13 @@ c           write (iout,*) "c(",j,i,") =",c(j,i)
           thetaref(i)=theta(i)
           phiref(i)=phi(i)
         enddo
-        nres=nres_temp
+c        nres=nres_temp
         if (waga_dist.ne.0.0d0) then
-          ii=0
-          do i = nnt,nct-2 
-            do j=i+2,nct 
+          ii=ii_old
+c          do i = nnt,nct-2 
+          do i = nnt_chain,nct_chain-2 
+c            do j=i+2,nct 
+            do j=i+2,nct_chain
 
               x12=c(1,i)-c(1,j)
               y12=c(2,i)-c(2,j)
@@ -3755,9 +3807,9 @@ c              write (iout,*) k,i,j,distal,dist2_cut
 
 c             write (iout,*) "k",k
 c             write (iout,*) "i",i," j",j," constr_homology",
-c    &                       constr_homology
-              ires_homo(ii)=i
-              jres_homo(ii)=j
+c     &                       constr_homology
+              ires_homo(ii)=i+chain_border1(1,igr)-1
+              jres_homo(ii)=j+chain_border1(1,igr)-1
               odl(k,ii)=distal
               if (read2sigma) then
                 sigma_odl(k,ii)=0
 c     Theta, dihedral and SC retraints
 c
         if (waga_angle.gt.0.0d0) then
-          do i = nnt+3,nct 
+          do i = nnt_chain+3,nct_chain
+            iii=i+chain_border1(1,igr)-1
             if (idomain(k,i).eq.0) then 
 c               sigma_dih(k,i)=0.0
                cycle
             endif
-            dih(k,i)=phiref(i)
-            sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
+            dih(k,iii)=phiref(i)
+            sigma_dih(k,iii)=
+     &         (rescore(k,i)+rescore(k,i-1)+
      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
 c            write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
 c     &       " sigma_dihed",sigma_dih(k,i)
-            if (sigma_dih(k,i).ne.0)
-     &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
+            if (sigma_dih(k,iii).ne.0)
+     &       sigma_dih(k,iii)=1.0d0/(sigma_dih(k,iii)*sigma_dih(k,iii))
           enddo
-          lim_dih=nct-nnt-2 
+c          lim_dih=nct-nnt-2 
         endif
 
         if (waga_theta.gt.0.0d0) then
-          do i = nnt+2,nct
+          do i = nnt_chain+2,nct_chain
+             iii=i+chain_border1(1,igr)-1
              if (idomain(k,i).eq.0) then  
 c              sigma_theta(k,i)=0.0
               cycle
              endif
-             thetatpl(k,i)=thetaref(i)
-             sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
+             thetatpl(k,iii)=thetaref(i)
+             sigma_theta(k,iii)=(rescore(k,i)+rescore(k,i-1)+
      &                        rescore(k,i-2))/3.0
-             if (sigma_theta(k,i).ne.0)
-     &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
+             if (sigma_theta(k,iii).ne.0)
+     &       sigma_theta(k,iii)=1.0d0/
+     &       (sigma_theta(k,iii)*sigma_theta(k,iii))
           enddo
         endif
 
         if (waga_d.gt.0.0d0) then
-          do i = nnt,nct
+          do i = nnt_chain,nct_chain
+             iii=i+chain_border1(1,igr)-1
                if (itype(i).eq.10) cycle 
                if (idomain(k,i).eq.0 ) then 
 c                  sigma_d(k,i)=0.0
                   cycle
                endif
-               xxtpl(k,i)=xxref(i)
-               yytpl(k,i)=yyref(i)
-               zztpl(k,i)=zzref(i)
-               sigma_d(k,i)=rescore(k,i)
-               if (sigma_d(k,i).ne.0)
-     &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
-               if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
+               xxtpl(k,iii)=xxref(i)
+               yytpl(k,iii)=yyref(i)
+               zztpl(k,iii)=zzref(i)
+               sigma_d(k,iii)=rescore(k,i)
+               if (sigma_d(k,iii).ne.0)
+     &          sigma_d(k,iii)=1.0d0/(sigma_d(k,iii)*sigma_d(k,iii))
+c               if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
           enddo
         endif
       enddo ! l
 c remove distance restraints not used in any model from the list
 c shift data in all arrays
 c
+c      write (iout,*) "ii_old",ii_old
       if (waga_dist.ne.0.0d0) then
-        ii=0
+#ifdef DEBUG
+       write (iout,*) "Distance restraints from templates"
+       do iii=1,lim_odl
+       write(iout,'(4i5,100(2f8.2,1x,l1,4x))') 
+     &  iii,ii_in_use(iii),ires_homo(iii),jres_homo(iii),
+     &  (odl(ki,iii),1.0d0/dsqrt(sigma_odl(ki,iii)),l_homo(ki,iii),
+     &  ki=1,constr_homology)
+       enddo
+#endif
+        ii=ii_old
         liiflag=.true.
-        do i=nnt,nct-2 
-         do j=i+2,nct 
+        lfirst=.true.
+        do i=nnt_chain,nct_chain-2
+         do j=i+2,nct_chain
           ii=ii+1
-          if (ii_in_use(ii).eq.0.and.liiflag) then
+c          if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
+c     &            .and. distal.le.dist2_cut ) then
+c          write (iout,*) "i",i," j",j," ii",ii," i_in_use",ii_in_use(ii)
+c          call flush(iout)
+          if (ii_in_use(ii).eq.0.and.liiflag.or.
+     &     ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
             liiflag=.false.
-            iistart=ii
+            i10=ii
+            if (lfirst) then
+              lfirst=.false.
+              iistart=ii
+            else
+              if(i10.eq.lim_odl) i10=i10+1
+              do ki=0,i10-i01-1
+               ires_homo(iistart+ki)=ires_homo(ki+i01)
+               jres_homo(iistart+ki)=jres_homo(ki+i01)
+               ii_in_use(iistart+ki)=ii_in_use(ki+i01)
+               do k=1,constr_homology
+                odl(k,iistart+ki)=odl(k,ki+i01)
+                sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
+                l_homo(k,iistart+ki)=l_homo(k,ki+i01)
+               enddo
+              enddo
+              iistart=iistart+i10-i01
+            endif
           endif
-          if (ii_in_use(ii).ne.0.and..not.liiflag.or.
-     &                   .not.liiflag.and.ii.eq.lim_odl) then
-             if (ii.eq.lim_odl) then
-              iishift=ii-iistart+1
-             else
-              iishift=ii-iistart
-             endif
+          if (ii_in_use(ii).ne.0.and..not.liiflag) then
+             i01=ii
              liiflag=.true.
-             do ki=iistart,lim_odl-iishift
-              ires_homo(ki)=ires_homo(ki+iishift)
-              jres_homo(ki)=jres_homo(ki+iishift)
-              ii_in_use(ki)=ii_in_use(ki+iishift)
-              do k=1,constr_homology
-               odl(k,ki)=odl(k,ki+iishift)
-               sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
-               l_homo(k,ki)=l_homo(k,ki+iishift)
-              enddo
-             enddo
-             ii=ii-iishift
-             lim_odl=lim_odl-iishift
           endif
          enddo
         enddo
+        lim_odl=iistart-1
       endif
 
+      lll=lim_odl-ii_old
+
+      do i=2,nequiv(igr) 
+
+        ichain=iequiv(i,igr)
+
+        do j=nnt_chain,nct_chain
+          jj=j+chain_border1(1,ichain)-chain_border1(1,iequiv(1,igr))
+          do k=1,constr_homology
+            dih(k,jj)=dih(k,j)
+            sigma_dih(k,jj)=sigma_dih(k,j)
+            thetatpl(k,jj)=thetatpl(k,j)
+            sigma_theta(k,jj)=sigma_theta(k,j)
+            xxtpl(k,jj)=xxtpl(k,j)
+            yytpl(k,jj)=yytpl(k,j)
+            zztpl(k,jj)=zztpl(k,j)
+            sigma_d(k,jj)=sigma_d(k,j)
+          enddo
+        enddo
+
+        jj=chain_border1(1,ichain)-chain_border1(1,iequiv(i-1,igr)) 
+c        write (iout,*) "igr",igr," i",i," ichain",ichain," jj",jj
+        do j=ii_old+1,lim_odl
+          ires_homo(j+lll)=ires_homo(j)+jj
+          jres_homo(j+lll)=jres_homo(j)+jj
+          do k=1,constr_homology
+            odl(k,j+lll)=odl(k,j)
+            sigma_odl(k,j+lll)=sigma_odl(k,j)
+            l_homo(k,j+lll)=l_homo(k,j)
+          enddo
+        enddo
+
+        ii_old=ii_old+lll
+        lim_odl=lim_odl+lll
+
+      enddo
+
+      ENDDO ! IGR
+
+      if (waga_angle.gt.0.0d0) lim_dih=nct-nnt-2 
+      nres=nres_temp
+      itype=itype_temp
+
       return
    10 stop "Error in fragment file"
       end
index 4b7b763..5573f05 100644 (file)
@@ -21,12 +21,12 @@ c      print *,'just initialize'
 c      print *,fail
       s1=0.0
       s2=0.0
-      print *,s1,s2
+c      print *,s1,s2
       do 1 i=1,3
-      print *, i2,i3,i4
+c      print *, i2,i3,i4
       zi=c(i,i2)-c(i,i3)
       ui=c(i,i4)-c(i,i3)
-      print *,zi,ui
+c      print *,zi,ui
       s1=s1+zi*zi
       s2=s2+ui*ui
       z(i)=zi
@@ -41,7 +41,7 @@ c      print *,fail
       write(iout,1000) i3,i4,i1
       fail=.true.
       return
-      print *,'two if pass'
+c      print *,'two if pass'
     4 s1=1.0/s1
       s2=1.0/s2
       v1=z(2)*u(3)-z(3)*u(2)
index c83e9ce..946fb58 100644 (file)
@@ -1,5 +1,8 @@
       subroutine friction_force
       implicit none
+#ifdef MPI
+      include 'mpif.h'
+#endif
       include 'DIMENSIONS'
       include 'COMMON.VAR'
       include 'COMMON.CHAIN'
       
       logical lprn /.false./, checkmode /.false./
 #ifdef FIVEDIAG
+#ifdef TIMING
+      include 'COMMON.TIME1'
+      double precision time01
+#endif
 c Here accelerations due to friction forces are computed right after forces.
-      d_t_work=0.0d0
+      d_t_work(:6*nres)=0.0d0
       do j=1,3
         v_work(j,1)=d_t(j,0)
         v_work(j,nnt)=d_t(j,0)
@@ -86,8 +93,14 @@ c          inct=chain_border(2,1)
           write (iout,'(f10.5)') (vvec(i),i=iposc,ind)
 #endif
 c          write (iout,*) "chain",i," ind",ind," n",n
+#ifdef TIMING
+          time01=MPI_Wtime()
+#endif
           call fivediagmult(n,DMfric(iposc),DU1fric(iposc),
      &     DU2fric(iposc),vvec(iposc),rs)
+#ifdef TIMING
+          time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01
+#endif
 #ifdef DEBUG
           write (iout,*) "rs"
           write (iout,'(f10.5)') (rs(i),i=1,n)
@@ -104,7 +117,7 @@ c     &       "index",3*(i-1)+j,"rs",rs(i-iposc+1)
       write (iout,'(3f10.5)') (fric_work(j),j=1,dimen3)
 #endif
 #else
-      do i=0,MAXRES2
+      do i=0,2*nres
         do j=1,3
           friction(j,i)=0.0d0
         enddo
@@ -277,7 +290,7 @@ c-----------------------------------------------------
       integer ichain,innt,inct,iposc
 #endif
 
-      do i=0,MAXRES2
+      do i=0,2*nres
         do j=1,3
           stochforc(j,i)=0.0d0
         enddo
@@ -567,9 +580,9 @@ C      gamsc(ntyp1)=1.0d0
         enddo
       endif
 #ifdef FIVEDIAG
-      DMfric=0.0d0
-      DU1fric=0.0d0
-      DU2fric=0.0d0
+      DMfric(:2*nres)=0.0d0
+      DU1fric(:2*nres)=0.0d0
+      DU2fric(:2*nres)=0.0d0
       ind=1
       do ichain=1,nchain
         innt=chain_border(1,ichain)
index 7bd51b8..98ce59a 100644 (file)
@@ -50,6 +50,18 @@ C Calculate the initial time, if it is not zero (e.g. for the SUN).
       time_fricmatmult=0.0d0
       time_fsample=0.0d0
       time_SAXS=0.0d0
+      time_list=0.0d0
+      time_evdw=0.0d0
+      time_evdw_short=0.0d0
+      time_evdw_long=0.0d0
+      time_eelec=0.0d0
+      time_eelec_short=0.0d0
+      time_eelec_long=0.0d0
+      time_escp=0.0d0
+      time_escp_short=0.0d0
+      time_escp_long=0.0d0
+      time_escpsetup=0.0d0
+      time_escpcalc=0.0d0
 #endif
 cd    print *,' in SET_TIMERS stime=',stime
       return 
@@ -287,6 +299,7 @@ C---------------------------------------------------------------------------
       include 'COMMON.IOUNITS'
       include 'COMMON.TIME1'
       include 'COMMON.SETUP'
+      include 'COMMON.MD'
 #ifdef MPI
       time1=MPI_WTIME()
          write (iout,'(80(1h=)/a/(80(1h=)))') 
@@ -318,6 +331,30 @@ C---------------------------------------------------------------------------
      &      time_bcast+time_reduce+time_gather+time_scatter+
      &      time_sendrecv+time_barrier_g+time_barrier_e+time_bcastc
          write (*,*) "Processor",fg_rank,myrank," enecalc",time_enecalc
+#ifdef TIMING_ENE
+         if (RESPA) then
+           write (*,*) "Processor",fg_rank,myrank," evdw_long",
+     &       time_evdw_long
+           write (*,*) "Processor",fg_rank,myrank," evdw_short",
+     &       time_evdw_short
+           write (*,*) "Processor",fg_rank,myrank," eelec_long",
+     &       time_eelec_long
+           write (*,*) "Processor",fg_rank,myrank," eelec_short",
+     &       time_eelec_short
+           write (*,*) "Processor",fg_rank,myrank," escp_long",
+     &       time_escp_long
+           write (*,*) "Processor",fg_rank,myrank," escp_short",
+     &       time_escp_short
+         else
+           write (*,*) "Processor",fg_rank,myrank," evdw",time_evdw
+           write (*,*) "Processor",fg_rank,myrank," eelec",time_eelec
+           write (*,*) "Processor",fg_rank,myrank," escp",time_escp
+           write (*,*) "Processor",fg_rank,myrank," escpsetup",
+     &      time_escpsetup
+           write (*,*) "Processor",fg_rank,myrank," escpcalc",
+     &      time_escpcalc
+         endif
+#endif
          write (*,*) "Processor",fg_rank,myrank," sumene",time_sumene
          write (*,*) "Processor",fg_rank,myrank," intfromcart",
      &     time_intfcart
index 978fa59..2e0ebaf 100644 (file)
@@ -213,11 +213,11 @@ c      include 'COMMON.CONTACTS'
       include 'COMMON.REMD'
       include 'COMMON.MD'
       include 'COMMON.SBRIDGE'
-      integer i,icall,iretcode,nfun
+      integer i,it,icall,iretcode,nfun
       common /srutu/ icall
       integer nharp,iharp(4,maxres/3)
       integer nft_sc
-      logical fail
+      logical fail,secondary_str /.true./
       double precision energy(0:n_ene),etot,etota
       double precision energy_long(0:n_ene),energy_short(0:n_ene)
       double precision rms,frac,frac_nn,co
@@ -243,6 +243,24 @@ c      include 'COMMON.CONTACTS'
       write (iout,*) "Energy evaluation/minimization"
       call chainbuild_cart
 c      print *,'dc',dc(1,0),dc(2,0),dc(3,0)
+      if (nran_start.gt.0) then
+        write (iout,*) 
+     & "Chains will be regenerated starting from residue",nran_start
+        do it=1,100
+        call gen_rand_conf_mchain(nran_start,*10)
+        write (iout,*) "Conformation successfully generated",it
+        goto 11
+   10   write (iout,*) "Problems with regenerating chains",it
+        enddo
+   11   continue
+        write (iout,*) "Cartesian coords after chain rebuild"
+        call cartprint
+        call chainbuild_cart
+        write (iout,*) "Cartesian coords after chainbuild_ecart"
+        call cartprint
+        call int_from_cart1(.false.)
+        call intout
+      endif
       if (split_ene) then
        print *,"Processor",myrank," after chainbuild"
        icall=1
@@ -274,9 +292,11 @@ c      print *,"after etotal"
       etota = energy(0)
       etot =etota
       call enerprint(energy(0))
+      if (secondary_str) then
       call hairpin(.true.,nharp,iharp)
 c        print *,'after hairpin'
       call secondary2(.true.)
+      endif
 c        print *,'after secondary'
       if (minim) then
 crc overlap test
@@ -292,6 +312,9 @@ crc overlap test
           write (iout,*) 'Calling OVERLAP_SC'
           call overlap_sc(fail)
           write (iout,*) "After overlap_sc"
+c        cartname=prefix(:ilen(prefix))//'.x'
+c        potE=etot
+c        call cartoutx(0.0d0)
         endif 
 
         if (searchsc) then 
@@ -330,10 +353,12 @@ c          print *,'Calling MINIMIZE.'
 #endif
         print *,'# eval/s',evals
         print *,'refstr=',refstr
+        if (secondary_str) then
         call hairpin(.false.,nharp,iharp)
-        print *,'after hairpin'
+c        print *,'after hairpin'
         call secondary2(.true.)
-        print *,'after secondary'
+c        print *,'after secondary'
+        endif
         call etotal(energy(0))
         etot = energy(0)
         call enerprint(energy(0))