Adam's wham update
authorCezary Czaplewski <czarek@chem.univ.gda.pl>
Tue, 16 Feb 2021 13:54:17 +0000 (14:54 +0100)
committerCezary Czaplewski <czarek@chem.univ.gda.pl>
Tue, 16 Feb 2021 13:54:17 +0000 (14:54 +0100)
19 files changed:
source/wham/src-HCD/COMMON.CHAIN
source/wham/src-HCD/COMMON.DFA
source/wham/src-HCD/COMMON.ENERGIES
source/wham/src-HCD/DIMENSIONS
source/wham/src-HCD/DIMENSIONS.FREE
source/wham/src-HCD/chain_symmetry.F
source/wham/src-HCD/cxread.F
source/wham/src-HCD/elecont.f
source/wham/src-HCD/enecalc1.F
source/wham/src-HCD/energy_p_new.F
source/wham/src-HCD/geomout.F
source/wham/src-HCD/include_unres/COMMON.SBRIDGE
source/wham/src-HCD/include_unres/COMMON.SETUP
source/wham/src-HCD/initialize_p.F
source/wham/src-HCD/make_ensemble1.F
source/wham/src-HCD/read_constr_homology.F
source/wham/src-HCD/secondary.f
source/wham/src-HCD/ssMD.F
source/wham/src-HCD/wham_calc1.F

index 7b79a58..6b453af 100644 (file)
@@ -3,7 +3,7 @@
      & tabpermchain,nchain ,npermchain,ireschain,iz_sc,nres_chomo
       double precision c,cref,crefjlee,dc,xloc,xrot,dc_norm,t,r,prod,rt,
      & rmssing,anatemp,chomo
-      common /chain/ c(3,maxres2+2),dc(3,maxres2),xloc(3,maxres),
+      common /chain/ c(3,maxres2+2),dc(3,0:maxres2),xloc(3,maxres),
      & xrot(3,maxres),dc_norm(3,maxres2),nres,nres0
       common /rotmat/ t(3,3,maxres),r(3,3,maxres),prod(3,3,maxres),
      & rt(3,3,maxres) 
index 064a7ce..782e8c4 100644 (file)
@@ -11,7 +11,7 @@ C Total : ~ 11 * Nres restraints
 C
 C
       INTEGER IDFAMAX,IDFAMX2,IDFACMD,IDMAXMIN, MAXN
-      PARAMETER(IDFAMAX=10000,IDFAMX2=1000,IDFACMD=500,IDMAXMIN=500)
+      PARAMETER(IDFAMAX=25000,IDFAMX2=1000,IDFACMD=500,IDMAXMIN=500)
       PARAMETER(MAXN=4)
       real*8 wwdist,wwangle,wwnei
       parameter(wwdist=1.0d0,wwangle=1.0d0,wwnei=1.0d0)
index 2d40a95..e0ed696 100644 (file)
@@ -1,4 +1,4 @@
       double precision potE(MaxStr_Proc,Max_Parm),entfac(MaxStr_Proc),
-     & q(MaxQ+2,MaxStr_Proc),enetb(max_ene,MaxStr_Proc,Max_Parm)
+     & q(MaxQ+6,MaxStr_Proc),enetb(max_ene,MaxStr_Proc,Max_Parm)
       integer einicheck
       common /energies/ potE,entfac,q,enetb,einicheck
index fb93081..49e50f5 100644 (file)
@@ -14,8 +14,8 @@ c      parameter (max_cg_procs=maxprocs)
 C Max. number of AA residues
       integer maxres
 c      parameter (maxres=250)
-c      parameter (maxres=1200)
-      parameter (maxres=10000)
+      parameter (maxres=1200)
+c      parameter (maxres=20000)
 C Max. number of cysteines and other bridging residues
       integer max_cyst
       parameter (max_cyst=100)
index 7a397d9..c42f2d7 100644 (file)
@@ -3,7 +3,7 @@
       integer MaxR,MaxT_h,maxHdim
       integer MaxSlice
       parameter (Max_Parm=5)
-      parameter (MaxQ=4,MaxQ1=MaxQ+2)
+      parameter (MaxQ=4,MaxQ1=MaxQ+6)
       parameter(MaxR=8,MaxT_h=36)
       parameter(MaxSlice=40)
       parameter(maxHdim=200)
index 1406d1d..8c36855 100644 (file)
@@ -7,6 +7,7 @@ c
       implicit none
       include "DIMENSIONS"
       include "COMMON.IOUNITS"
+      include "COMMON.CONTROL"
       integer nchain,nres,itype(nres),chain_border(2,maxchain),
      &  chain_length(nchain),itemp(maxchain),
      &  npermchain,tabpermchain(maxchain,maxperm),
@@ -42,6 +43,7 @@ c
         nchain_group=nchain_group+1
         iieq=1
         iequiv(iieq,nchain_group)=i
+        if (symetr.eq.1) then
         do j=i+1,nchain 
           if (iflag(j).gt.0.or.chain_length(i).ne.chain_length(j)) cycle
 c          k=0
@@ -57,6 +59,7 @@ c            k=k+1
           iieq=iieq+1
           iequiv(iieq,nchain_group)=j
         enddo
+        endif
         nequiv(nchain_group)=iieq
       enddo
       write(iout,*) "Number of equivalent chain groups:",nchain_group
index 46867c5..2cfb938 100644 (file)
@@ -174,7 +174,7 @@ c      call flush(iout)
 c      write (iout,*) "Before boxshift"
 c      call flush(iout)
 c Box shift
-      call oligomer
+c      call oligomer
 c      write (iout,*) "After oligomer"
 c      call flush(iout)
       do i=1,nres
index fb105a4..86db2df 100644 (file)
      &  eesij,ees,evdw,ene, rij,zj_temp,xj_temp,yj_temp,
      & sscale,sscagrad,dist_temp,xj_safe,yj_safe,zj_safe,dist_init
       double precision elpp6c(2,2),elpp3c(2,2),ael6c(2,2),ael3c(2,2),
-     &  appc(2,2),bppc(2,2)
+     &  appc(2,2),bppc(2,2),epp_(2,2),rpp_(2,2)
       double precision elcutoff,elecutoff_14
       integer ncont,icont(2,maxcont),xshift,yshift,zshift,isubchap
       double precision econt(maxcont)
+      double precision boxshift
 *
 * Load the constants of peptide bond - peptide bond interactions.
 * Type 1 - ordinary peptide bond, type 2 - alkylated peptide bond (e.g.
@@ -29,8 +30,8 @@
 *
 * as of 7/06/91.
 *
-c      data epp    / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/
-c      data rpp    / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
+c      data epp_   / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/
+      data rpp_   / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
       data elpp6c  /-0.2379d0,-0.2056d0,-0.2056d0,-0.0610d0/
       data elpp3c  / 0.0503d0, 0.0000d0, 0.0000d0, 0.0692d0/
       data elcutoff /-0.3d0/,elecutoff_14 /-0.5d0/
@@ -40,7 +41,7 @@ c      data rpp    / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
      &  "Constants of electrostatic interaction energy expression."
       do i=1,2
         do j=1,2
-        rri=rpp(i,j)**6
+        rri=rpp_(i,j)**6
         appc(i,j)=epp(i,j)*rri*rri 
         bppc(i,j)=-2.0*epp(i,j)*rri
         ael6c(i,j)=elpp6c(i,j)*4.2**6
@@ -62,12 +63,8 @@ c      data rpp    / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
         xmedi=xi+0.5*dxi
         ymedi=yi+0.5*dyi
         zmedi=zi+0.5*dzi
-          xmedi=mod(xmedi,boxxsize)
-          if (xmedi.lt.0) xmedi=xmedi+boxxsize
-          ymedi=mod(ymedi,boxysize)
-          if (ymedi.lt.0) ymedi=ymedi+boxysize
-          zmedi=mod(zmedi,boxzsize)
-          if (zmedi.lt.0) zmedi=zmedi+boxzsize
+        call to_box(xmedi,ymedi,zmedi)
+c        write (iout,*) "i",xmedi,ymedi,zmedi
         do 4 j=i+2,ien-1
           jj=iperm(j,ipermmin)
           ind=ind+1
@@ -86,46 +83,13 @@ c      data rpp    / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
           xj=c(1,jj)+0.5*dxj
           yj=c(2,jj)+0.5*dyj
           zj=c(3,jj)+0.5*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
-      isubchap=0
-      do xshift=-1,1
-      do yshift=-1,1
-      do zshift=-1,1
-          xj=xj_safe+xshift*boxxsize
-          yj=yj_safe+yshift*boxysize
-          zj=zj_safe+zshift*boxzsize
-          dist_temp=(xj-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
-            isubchap=1
-          endif
-       enddo
-       enddo
-       enddo
-       if (isubchap.eq.1) then
-          xj=xj_temp-xmedi
-          yj=yj_temp-ymedi
-          zj=zj_temp-zmedi
-       else
-          xj=xj_safe-xmedi
-          yj=yj_safe-ymedi
-          zj=zj_safe-zmedi
-       endif
+c          write (iout,*) "j",xj,yj,zj
+          call to_box(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
-            sss=sscale(sqrt(rij))
-            sssgrad=sscagrad(sqrt(rij))
           rrmij=1.0/(xj*xj+yj*yj+zj*zj)
           rmij=sqrt(rrmij)
           r3ij=rrmij*rmij
@@ -152,6 +116,7 @@ c      data rpp    / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
           endif
           ees=ees+eesij
           evdw=evdw+evdwij*sss
+c          write (iout,*) "i",i," j",j," rij",dsqrt(rij)," eesij",eesij
     4   continue
     1 continue
       if (lprint) then
index b64de48..2edf349 100644 (file)
@@ -46,13 +46,31 @@ c      double precision tole /1.0d-1/
       double precision tt
       integer snk_p(MaxR,MaxT_h,Max_parm)
       logical lerr
+      integer ncont,icont(2,maxcont),isecstr(maxres)
       character*256 bprotfile_temp
+      double precision totlength
       call opentmp(islice,ientout,bprotfile_temp)
       iii=0
       ii=0
       errmsg_count=0
 c      write (iout,*) "enecalc: nparmset ",nparmset
 c      write (iout,*) "enecalc: tormode ",tor_mode
+      write (iout,*) "ns",ns," dyn_ss",dyn_ss,(iss(i),i=1,ns)
+      if (ns.gt.0.and.dyn_ss) then
+          do i=nss+1,nhpb
+            ihpb(i-nss)=ihpb(i)
+            jhpb(i-nss)=jhpb(i)
+            forcon(i-nss)=forcon(i)
+            dhpb(i-nss)=dhpb(i)
+          enddo
+          nhpb=nhpb-nss
+          nss=0
+          call hpb_partition
+          do i=1,ns
+            dyn_ss_mask(iss(i))=.true.
+          enddo
+      endif
+      write (iout,*) "dyn_ss_mask",(dyn_ss_mask(i),i=1,nres)
 #ifdef MPI
       do iparm=1,nParmSet
         do ib=1,nT_h(iparm)
@@ -94,7 +112,7 @@ c      write (iout,*) "enecalc: tormode ",tor_mode
            anatemp= 1.0d0/(beta_h(ib,ipar)*1.987D-3)
            q(nQ+1,iii+1)=rmsnat(iii+1,ipermin)
          endif
-         q(nQ+2,iii+1)=gyrate(iii+1)
+c         write (iout,*) iii+1,q(nQ+3,iii+1),q(nQ+4,iii+1),q(nQ+5,iii+1)
 c        fT=T0*beta_h(ib,ipar)*1.987D-3
 c        ft=2.0d0/(1.0d0+1.0d0/(T0*beta_h(ib,ipar)*1.987D-3))
         if (rescale_mode.eq.1) then
@@ -158,9 +176,9 @@ c     &   " kfac",kfac,"quot",quot," fT",fT
      &      wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
      &      wtor_d,wsccor,wbond
 #endif
-C        write (iout,*) "tuz przed energia"
+c        write (iout,*) "tuz przed energia"
         call etotal(energia(0),fT)
-C        write (iout,*) "tuz za energia"
+c        write (iout,*) "tuz za energia"
 #ifdef DEBUG
         write (iout,*) "Conformation",i
         write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres),
@@ -215,8 +233,8 @@ c     &       eini-energia(0)
           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres),
      &                            ((c(l,k+nres),l=1,3),k=nnt,nct)
 c              call intout
-              call pdbout(indstart(me1)+iii,
-     & 1.0d0/(1.987D-3*beta_h(ib,ipar)),energia(0),eini,0.0d0,0.0d0)
+c              call pdbout(indstart(me1)+iii,
+c     & 1.0d0/(1.987D-3*beta_h(ib,ipar)),energia(0),eini,0.0d0,0.0d0)
               call enerprint(energia(0),fT)
               errmsg_count=errmsg_count+1
               if (errmsg_count.gt.maxerrmsg_count) 
@@ -260,6 +278,30 @@ c          call enerprint(energia(0),fT)
         endif
 
         enddo ! iparm
+        q(nQ+2,iii+1)=gyrate(iii+1)
+c 8/28/2020 Adam - determine the fraction of secondary structures.
+        call elecont(.false.,ncont,icont,nnt,nct-1,1)
+        call secondary2(.false.,.false.,ncont,icont,isecstr)
+#ifdef DEBUG
+        write (iout,*) "secondary structure"
+        write (iout,'(80i1)') (isecstr(k),k=1,nres)
+#endif
+        q(nQ+3,iii+1)=0.0d0
+        q(nQ+4,iii+1)=0.0d0
+        q(nQ+5,iii+1)=0.0d0
+        totlength=0.0d0
+        do k=nnt,nct
+          if (itype(k).eq.ntyp1) cycle
+          totlength=totlength+1.0d0
+          l=isecstr(k)
+          q(nQ+3+l,iii+1)=q(nQ+3+l,iii+1)+1.0d0
+        enddo
+        q(nQ+3,iii+1)=q(nQ+3,iii+1)/totlength
+        q(nQ+4,iii+1)=q(nQ+4,iii+1)/totlength
+        q(nQ+5,iii+1)=q(nQ+5,iii+1)/totlength
+c        write (iout,*) "iii",iii," nssbond",nssbond,nss
+c        q(nQ+6,iii+1)=nssbond
+        q(nQ+6,iii+1)=nss
 
         iii=iii+1
         if (q(1,iii).le.0.0d0 .and. indpdb.gt.0)
index e72d558..3a83918 100644 (file)
@@ -56,7 +56,7 @@ C
        call set_shield_fac2
       endif
       call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
-C            write(iout,*) 'po eelec'
+c            write(iout,*) 'po eelec eello_turn4',eello_turn4
 
 C Calculate excluded-volume interaction energy between peptide groups
 C and side chains.
@@ -2273,6 +2273,7 @@ c        write(iout,*) "JESTEM W PETLI"
         call eelecij(i,i+3,ees,evdw1,eel_loc)
         if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1) 
      &   call eturn4(i,eello_turn4)
+c        write (iout,*) "i",i," eello_turn4",eello_turn4
 #ifdef FOURBODY
         num_cont_hb(i)=num_conti
 #endif
@@ -3581,6 +3582,8 @@ C Third- and fourth-order contributions from turns
       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
      &    dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
      &    num_conti,j1,j2
+      double precision sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij
+      common /lipcalc/ sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij
       j=i+3
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 C
index 097040f..ed33cc7 100644 (file)
@@ -9,7 +9,9 @@
       include 'COMMON.HEADER'
       include 'COMMON.SBRIDGE'
       character*50 tytul
-      character*1 chainid(10) /'A','B','C','D','E','F','G','H','I','J'/
+      character*1 chainid(32) /'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',
+     & '1','2','3','4','5','6'/
       dimension ica(maxres)
       write(ipdb,'("REMARK CONF",i8," TEMPERATURE",f7.1," RMS",0pf7.2)') 
      &  ii,temp,rmsdev
       iatom=0
       ichain=1
       ires=0
+      iti_prev=0
       do i=nnt,nct
         iti=itype(i)
         if (iti.eq.ntyp1) then
-          ichain=ichain+1
           ires=0
-          write (ipdb,'(a)') 'TER'
+          if (iti_prev.ne.ntyp1) then
+            write (ipdb,'(a)') 'TER'
+            ichain=ichain+1
+          endif
         else
         ires=ires+1
         iatom=iatom+1
         ica(i)=iatom
         write (ipdb,10) iatom,restyp(iti),chainid(ichain),
-     &     ires,(c(j,i),j=1,3)
+     &     ires,(c(j,i),j=1,3),1.0d0
         if (iti.ne.10) then
           iatom=iatom+1
           write (ipdb,20) iatom,restyp(iti),chainid(ichain),
-     &      ires,(c(j,nres+i),j=1,3)
+     &      ires,(c(j,nres+i),j=1,3),1.0d0
         endif
         endif
+        iti_prev=iti
       enddo
       write (ipdb,'(a)') 'TER'
       do i=nnt,nct-1
         write (ipdb,30) ica(nct),ica(nct)+1
       endif
       do i=1,nss
-        if (dyn_ss) then
-         write (iunit,30) ica(idssb(i))+1,ica(jdssb(i))+1
-        else
-         write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
-        endif
+        write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
       enddo
       write (ipdb,'(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)
+  10  FORMAT ('ATOM',I7,'  CA  ',A3,1X,A1,I4,4X,3F8.3,2f6.2)
+  20  FORMAT ('ATOM',I7,'  CB  ',A3,1X,A1,I4,4X,3F8.3,2f6.2)
   30  FORMAT ('CONECT',8I5)
       return
       end
index a313d8f..d2a41e1 100644 (file)
@@ -24,8 +24,9 @@
      & link_end_peak
       double precision Ht,dyn_ssbond_ij,dtriss,atriss,btriss,ctriss
       logical dyn_ss,dyn_ss_mask
+      integer nssbond
       common /dyn_ssbond/ dtriss,atriss,btriss,ctriss,Ht,
      &  dyn_ssbond_ij(max_cyst,max_cyst),
-     &  idssb(maxss),jdssb(maxss)
+     &  idssb(maxss),jdssb(maxss),nssbond
       common /dyn_ss_logic/
      &  dyn_ss,dyn_ss_mask(maxres)
index 5039116..1edc7c3 100644 (file)
@@ -1,14 +1,14 @@
       integer king,idint,idreal,idchar,is_done
       parameter (king=0,idint=1105,idreal=1729,idchar=1597,is_done=1)
       integer me,cg_rank,fg_rank,fg_rank1,nodes,Nprocs,nfgtasks,kolor,
-     & koniec(0:maxprocs-1),WhatsUp,ifinish(maxprocs-1),CG_COMM,FG_COMM,
+     & koniec(0:maxprocs-1),ifinish(maxprocs-1),CG_COMM,FG_COMM,
      & FG_COMM1,CONT_FROM_COMM,CONT_TO_COMM,lentyp(0:maxprocs-1),
      & kolor1,key1,nfgtasks1,MyRank,
      & max_gs_size
       logical yourjob, finished, cgdone
       common/setup/me,MyRank,cg_rank,fg_rank,fg_rank1,nodes,Nprocs,
      & nfgtasks,nfgtasks1,
-     & max_gs_size,kolor,koniec,WhatsUp,ifinish,CG_COMM,FG_COMM,
+     & max_gs_size,kolor,koniec,ifinish,CG_COMM,FG_COMM,
      & FG_COMM1,CONT_FROM_COMM,CONT_TO_COMM,lentyp
       integer MPI_UYZ,MPI_UYZGRAD,MPI_MU,MPI_MAT1,MPI_MAT2,
      & MPI_THET,MPI_GAM,
index 8217663..0cf093a 100644 (file)
@@ -22,6 +22,7 @@ C
       include "COMMON.NAMES"
       include "COMMON.TIME1"
       include "COMMON.TORCNSTR"
+      include "COMMON.SETUP"
 C
 C The following is just to define auxiliary variables used in angle conversion
 C
@@ -237,6 +238,7 @@ C Set timers and counters for the respective routines
       n_viol = 0
       n_gviol = 0
       n_map = 0
+      nfgtasks = 1
 #ifndef SPLITELE
       nprint_ene=nprint_ene-1
 #endif
index a07dbeb..f18dd2b 100644 (file)
@@ -371,7 +371,7 @@ c          write (iout,*) "qfree",qfree
      &        ctemper(:ilen(ctemper))//"pdb" 
           endif
           open(ipdb,file=pdbname)
-          write (iout,*) "Before reading nlist",nlist
+c          write (iout,*) "Before reading nlist",nlist
           do i=1,nlist
             read (ientout,rec=iperm(i))
      &        ((csingle(l,k),l=1,3),k=1,nres),
index fa81b80..03d7968 100644 (file)
@@ -259,8 +259,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
index 4088831..d2e9271 100644 (file)
@@ -401,8 +401,8 @@ cd      write (iout,*) '------- looking for parallel beta -----------'
       do i=1,ncont
         i1=icont(1,i)
         j1=icont(2,i)
-        if (i1.ge.nstart_sup .and. i1.le.nend_sup 
-     &     .and. j1.gt.nstart_sup .and. j1.le.nend_sup) then
+c        if (i1.ge.nstart_sup .and. i1.le.nend_sup 
+c     &     .and. j1.gt.nstart_sup .and. j1.le.nend_sup) then
 cd        write (iout,*) "parallel",i1,j1
         if(j1-i1.gt.5 .and. freeres(i1,j1,nsec,isec)) then
           ii1=i1
@@ -469,7 +469,7 @@ cd            write (iout,*) i1,j1,not_done
      &          "SetNeigh",ii1-1,i1-1,jj1-1,j1-1
            endif
           endif
-        endif
+c        endif
         endif ! i1.ge.nstart_sup .and. i1.le.nend_sup .and. i2.gt.nstart_sup .and. i2.le.nend_sup
       enddo
 
index 4ce1b3d..d9b9df7 100644 (file)
@@ -144,13 +144,13 @@ c-------TESTING CODE
       double precision echeck(-1:1),deps,ssx0,ljx0,xi,yi,zi
 c-------END TESTING CODE
 
-
+      nssbond=0
       i=resi
       j=resj
       ici=icys(i)
       icj=icys(j)
       if (ici.eq.0 .or. icj.eq.0) then
-        write (*,'(a,i5,2a,a3,i5,5h and ,a3,i5)') 
+        write (iout,'(a,i5,2a,a3,i5,5h and ,a3,i5)') 
      &  "Attempt to create",
      &  " a disulfide link between non-cysteine residues ",restyp(i),i,
      &  restyp(j),j
@@ -276,6 +276,8 @@ c     Stop and plot energy and derivative as a function of distance
      &       ((ljxm-ssxm+2.0d0)*icheck)/nicheck+jcheck*deps
 c-------END TESTING CODE
 
+c      write (iout,'(2(a,i5),4(a,f7.2))') "resi",resi," resj",resj,
+c     &  " ljxm",ljxm," ljxs",ljxs," ssxm",ssxm," rij",rij
       if (rij.gt.ljxm) then
         havebond=.false.
         ljd=rij-ljXs
@@ -298,6 +300,8 @@ c-------END TESTING CODE
      &       -2.0D0*alf12*eps3der+sigder*sigsq_om12
       else if (rij.lt.ssxm) then
         havebond=.true.
+        nssbond=nssbond+1
+c        write (iout,*) "ssMD: nssbond",nssbond
         ssd=rij-ssXs
         eij=ssA*ssd*ssd+ssB*ssd+ssC
         eij=eij*sss        
@@ -309,6 +313,7 @@ c-------END TESTING CODE
         eom2= 2*akth*deltat2+pom1-om1*pom2
         eom12=pom2
       else
+c          nssbond=nssbond+1
         omega=v1ss+2.0d0*v2ss*cosphi+3.0d0*v3ss*cosphi*cosphi
 
         d_ssxm(1)=0.5D0*akct/ssA
@@ -497,8 +502,8 @@ cgrad        enddo
 cgrad      enddo
 
       do l=1,3
-        gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(k)
-        gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(k)
+        gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(l)
+        gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l)
       enddo
 
       return
@@ -556,7 +561,6 @@ c     Includes
       include 'COMMON.SBRIDGE'
       include 'COMMON.CHAIN'
       include 'COMMON.IOUNITS'
-C      include 'COMMON.SETUP'
 #ifndef CLUST
 #ifndef WHAM
 C      include 'COMMON.MD'
@@ -572,6 +576,7 @@ c     Local variables
       logical found
       integer i_newnss(1024),displ(0:1024)
       integer g_newihpb(maxdim_cont),g_newjhpb(maxdim_cont),g_newnss
+      nfgtasks=1
 
       allnss=0
       do i=1,ns-1
index 7e4512d..2d1d661 100644 (file)
@@ -639,7 +639,7 @@ c     &   " WHAM_COMM",WHAM_COMM
           sumE_p(i,iparm)=0.0d0
           sumEbis_p(i,iparm)=0.0d0
           sumEsq_p(i,iparm)=0.0d0
-          do j=1,nQ+2
+          do j=1,nQ+6
             sumQ_p(j,i,iparm)=0.0d0
             sumQsq_p(j,i,iparm)=0.0d0
             sumEQ_p(j,i,iparm)=0.0d0
@@ -654,7 +654,7 @@ c     &   " WHAM_COMM",WHAM_COMM
           sumE(i,iparm)=0.0d0
           sumEbis(i,iparm)=0.0d0
           sumEsq(i,iparm)=0.0d0
-          do j=1,nQ+2
+          do j=1,nQ+6
             sumQ(j,i,iparm)=0.0d0
             sumQsq(j,i,iparm)=0.0d0
             sumEQ(j,i,iparm)=0.0d0
@@ -854,7 +854,7 @@ c          call restore_parm(iparm)
             sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
             sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
             sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
-            do j=1,nQ+2
+            do j=1,nQ+6
               sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
               sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
               sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm)
@@ -865,7 +865,7 @@ c          call restore_parm(iparm)
             sumE(k,iparm)=sumE(k,iparm)+etot*weight
             sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
             sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
-            do j=1,nQ+2
+            do j=1,nQ+6
               sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
               sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
               sumEQ(j,k,iparm)=sumEQ(j,k,iparm)
@@ -1085,7 +1085,7 @@ c          call restore_parm(iparm)
      &    sumW(i,iparm)
         sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm)
      &    -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2)
-        do j=1,nQ+2
+        do j=1,nQ+6
           sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm)
           sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm)
      &     -sumQ(j,i,iparm)**2
@@ -1096,15 +1096,15 @@ c          call restore_parm(iparm)
      &   (startGridT+i*delta_T))+potEmin
         write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
      &   sumW(i,iparm),sumE(i,iparm)
-        write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
+        write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+6)
         write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
-     &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
+     &   (sumQsq(j,i,iparm),j=1,nQ+6),(sumEQ(j,i,iparm),j=1,nQ+6)
         write (iout,*) 
         write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
      &   sumW(i,iparm),sumE(i,iparm)
-        write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
+        write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+6)
         write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
-     &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
+     &   (sumQsq(j,i,iparm),j=1,nQ+6),(sumEQ(j,i,iparm),j=1,nQ+6)
         write (34,*) 
         call flush(34)
       enddo