changes
authorCezary Czaplewski <czarek@chem.univ.gda.pl>
Sun, 19 Jul 2020 17:32:07 +0000 (19:32 +0200)
committerCezary Czaplewski <czarek@chem.univ.gda.pl>
Sun, 19 Jul 2020 17:32:07 +0000 (19:32 +0200)
14 files changed:
source/unres/src-HCD-5D/COMMON.DFA
source/unres/src-HCD-5D/DIMENSIONS
source/unres/src-HCD-5D/boxshift.f
source/unres/src-HCD-5D/chain_symmetry.F
source/unres/src-HCD-5D/dfa.F
source/unres/src-HCD-5D/energy_p_new_barrier.F
source/unres/src-HCD-5D/readpdb-mult.F
source/unres/src-HCD-5D/readrtns_CSA.F
source/unres/src-HCD-5D/ssMD.F
source/wham/src-HCD/bxread.F
source/wham/src-HCD/cxread.F
source/wham/src-HCD/enecalc1.F
source/wham/src-HCD/readrtns.F
source/wham/src-HCD/xread.F

index 818c024..51a7af7 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 ed21dfe..4236776 100644 (file)
@@ -16,7 +16,7 @@ C Max. number of coarse-grain processors
       parameter (max_cg_procs=maxprocs)
 C Max. number of AA residues
       integer maxres
-      parameter (maxres=10000)
+      parameter (maxres=20000)
 C Max. number of AA residues per chain
       integer maxres_chain
       parameter (maxres_chain=1200)
@@ -53,14 +53,15 @@ C Max. number of contacts per residue
 c      parameter (maxconts=50)
 C Max. number of interactions within cutoff per residue
       integer maxint_res
-      parameter (maxint_res=200)
+      parameter (maxint_res=250)
 C Max. number od residues within distance cufoff from a given residue to 
 C     include in template-based/contact distance restraints.
       integer maxcont_res
       parameter (maxcont_res=200)
 C Max. number of distance/contact-distance restraints
       integer maxdim_cont
-      parameter (maxdim_cont=maxres*maxcont_res)
+c      parameter (maxdim_cont=maxres*maxcont_res)
+      parameter (maxdim_cont=maxres*1000)
 C Number of AA types (at present only natural AA's will be handled
       integer ntyp,ntyp1
       parameter (ntyp=24,ntyp1=ntyp+1)
index 29d3406..070d9c9 100644 (file)
@@ -72,11 +72,17 @@ c--------------------------------------------------------------------------
       subroutine lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
       implicit none
       include 'DIMENSIONS'
+      include 'COMMON.IOUNITS'
       include 'COMMON.CHAIN'
       double precision xi,yi,zi,sslipi,ssgradlipi
       double precision fracinbuf
       double precision sscalelip,sscagradlip
-
+#define DEBUG
+#ifdef DEBUG
+      write (iout,*) "bordlipbot",bordlipbot," bordliptop",bordliptop
+      write (iout,*) "buflipbot",buflipbot," lipbufthick",lipbufthick
+      write (iout,*) "xi yi zi",xi,yi,zi
+#endif
       if ((zi.gt.bordlipbot).and.(zi.lt.bordliptop)) then
 C the energy transfer exist
         if (zi.lt.buflipbot) then
@@ -97,5 +103,9 @@ C lipbufthick is thickenes of lipid buffore
         sslipi=0.0d0
         ssgradlipi=0.0
       endif
+#ifdef DEBUG
+      write (iout,*) "sslipi",sslipi," ssgradlipi",ssgradlipi
+#endif
+#undef DEBUG
       return
       end
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 3982b6d..1af0b44 100644 (file)
@@ -368,6 +368,7 @@ 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 6d5b25f..6d6a817 100644 (file)
@@ -2044,7 +2044,8 @@ c              write(iout,*) "PO ZWYKLE", evdwij
               if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') 
      &                        'evdw',i,j,evdwij,' ss'
 C triple bond artifac removal
-              do k=j+1,iend(i,iint) 
+c              do k=j+1,iend(i,iint) 
+              do k=j+1,nct
 C search over all next residues
                 if (dyn_ss_mask(k)) then
 C check if they are cysteins
index 76cb6b6..41fe7f6 100644 (file)
@@ -64,9 +64,9 @@ c        call flush(iout)
           iterter(ires_old-1)=1
           itype(ires_old)=ntyp1
           iterter(ires_old)=1
-c          ishift1=ishift1+1
+          ishift1=ishift1+1
           ibeg=2
-          write (iout,*) "Chain ended",ires,ishift,ires_old,ibeg
+          write (iout,*) "Chain ended",ires,ishift,ires_old
           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,7 +115,6 @@ 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
@@ -134,13 +133,10 @@ c              write (iout,*) "BEG ires",ires
             else if (ibeg.eq.2) then
 ! Start a new chain
               ishift=-ires_old+ires-1 !!!!!
-c              ishift1=ishift1-1    !!!!!
-c              write (iout,*) "New chain started",ires,ires_old,ishift,
-c     &           ishift1
+              ishift1=ishift1-1    !!!!!
+c              write (iout,*) "New chain started",ires,ishift,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)
@@ -163,8 +159,7 @@ 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)
-c            write (iout,*) "backbone ",atom
-c            write (iout,*) ires,res,(c(j,ires),j=1,3)
+!            write (iout,*) "backbone ",atom
 #ifdef DEBUG
             write (iout,'(i6,i3,2x,a,3f8.3)') 
      &      ires,itype(ires),res,(c(j,ires),j=1,3)
@@ -234,7 +229,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*(-e1(j)+e2(j))/sqrt(2.0d0)
+              c(j,i)=c(j,i+1)-1.9d0*e2(j)
             enddo
           else !unres_pdb
            do j=1,3
@@ -417,9 +412,8 @@ C and convert the peptide geometry into virtual-chain geometry.
       character*3 seq,res
       character*5 atom
       character*80 card
-      double precision sccor(3,50)
+      double precision sccor(3,20)
       integer rescode,iterter(maxres)
-      logical zero
       do i=1,maxres
          iterter(i)=0
       enddo
@@ -587,7 +581,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*(-e1(j)+e2(j))/sqrt(2.0d0)
+              c(j,i)=c(j,i+1)-1.9d0*e2(j)
             enddo
           else !unres_pdb
            do j=1,3
@@ -622,7 +616,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*(-e1(j)+e2(j))/sqrt(2.0d0)
+            c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
           enddo
         else
         do j=1,3
@@ -654,7 +648,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*(e1(j)-e2(j))/dsqrt(2.0d0)
+            c(j,1)=c(j,2)-1.9d0*e2(j)
           enddo
         else
         do j=1,3
@@ -682,18 +676,6 @@ 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.)
@@ -701,7 +683,6 @@ 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 eeaf74c..4fbc0f1 100644 (file)
@@ -741,7 +741,7 @@ C
       integer ilen
       external ilen
       integer iperm,tperm
-      integer i,j,ii,k,l,itrial,itmp,i1,i2,it1,it2,nres_temp
+      integer i,j,ii,k,l,itrial,itmp,i1,i2,it1,it2,nres_temp,itemp
       double precision sumv
 C
 C Read PDB structure if applicable
@@ -1184,8 +1184,7 @@ c      write (iout,*) "After read_dist_constr nhpb",nhpb
  335        continue
             unres_pdb=.false.
             nres_temp=nres
-c            call readpdb
-            call readpdb_template(nmodel_start+1)
+            call readpdb
             close(ipdbin)
             if (nres.ge.nres_temp) then
               nmodel_start=nmodel_start+1
@@ -1196,11 +1195,22 @@ c            call readpdb
                 enddo
               enddo
             else
-              if (me.eq.king .or. .not. out1file)
+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)') 
      &           "Different number of residues",nres_temp,nres,
      &           " model skipped."
             endif
+  116       continue
             nres=nres_temp
           enddo
   332     continue
index 26807a0..67e9f10 100644 (file)
@@ -100,6 +100,7 @@ C-----------------------------------------------------------------------------
       include 'COMMON.IOUNITS'
       include 'COMMON.CALC'
       include 'COMMON.NAMES'
+      include 'COMMON.SPLITELE'
 #ifndef CLUST
 #ifndef WHAM
       include 'COMMON.MD'
@@ -153,6 +154,8 @@ c-------END TESTING CODE
       j=resj
       ici=icys(i)
       icj=icys(j)
+c      write (iout,*) "dyn_ssbond",resi,resj,ici,icj
+c      call flush(iout)
       if (ici.eq.0 .or. icj.eq.0) then
 #ifdef MPI
         write (*,'(a,i5,2a,a3,i5,5h and ,a3,i5)') 
@@ -177,6 +180,8 @@ c-------END TESTING CODE
       yi=c(2,nres+i)
       zi=c(3,nres+i)
       call to_box(xi,yi,zi)
+c      write (iout,*) "After to_box i",xi,yi,zi
+c      call flush(iout)
 C define scaling factor for lipids
 
 C        if (positi.le.0) positi=positi+boxzsize
@@ -184,12 +189,18 @@ C        print *,i
 C first for peptide groups
 c for each residue check if it is in lipid or lipid water border area
       call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
+c      write (iout,*) "After lipid_layer"
+c      call flush(iout)
       itypj=itype(j)
       xj=c(1,nres+j)
       yj=c(2,nres+j)
       zj=c(3,nres+j)
       call to_box(xj,yj,zj)
+c      write (iout,*) "After to_box j",xj,yj,zj
+c      call flush(iout)
       call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+c      write (iout,*) "After lipid_layer"
+c      call flush(iout)
       aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
      &  +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
       bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
@@ -197,6 +208,8 @@ c for each residue check if it is in lipid or lipid water border area
       xj=boxshift(xj-xi,boxxsize)
       yj=boxshift(yj-yi,boxysize)
       zj=boxshift(zj-zi,boxzsize)
+c      write (iout,*) "After boxshift"
+c      call flush(iout)
       dxj=dc_norm(1,nres+j)
       dyj=dc_norm(2,nres+j)
       dzj=dc_norm(3,nres+j)
@@ -214,8 +227,8 @@ c for each residue check if it is in lipid or lipid water border area
 
       rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
       rij=dsqrt(rrij)  ! sc_angular needs rij to really be the inverse
-            sss=sscale((1.0d0/rij)/sigma(itypi,itypj))
-            sssgrad=sscagrad((1.0d0/rij)/sigma(itypi,itypj))
+            sss=sscale((1.0d0/rij)/sigma(itypi,itypj),r_cut_int)
+            sssgrad=sscagrad((1.0d0/rij)/sigma(itypi,itypj),r_cut_int)
 c     The following are set in sc_angular
 c      erij(1)=xj*rij
 c      erij(2)=yj*rij
@@ -223,7 +236,11 @@ c      erij(3)=zj*rij
 c      om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
 c      om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
 c      om12=dxi*dxj+dyi*dyj+dzi*dzj
+c      write (iout,*) "Calling sc_angular"
+c      call flush(iout)
       call sc_angular
+c      write (iout,*) "After sc_angular"
+c      call flush(iout)
       rij=1.0D0/rij  ! Reset this so it makes sense
 
       sig0ij=sigma(itypi,itypj)
@@ -265,6 +282,8 @@ c-------END TESTING CODE
 
 c-------TESTING CODE
 c     Stop and plot energy and derivative as a function of distance
+c      write (iout,*) "checkstop",checkstop
+c      call flush(iout)
       if (checkstop) then
         ssm=ssC-0.25D0*ssB*ssB/ssA
         ljm=-0.25D0*ljB*bb/aa
@@ -447,6 +466,8 @@ c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
 
       endif
 
+c      write (iout,*) "havebond",havebond
+c      call flush(iout)
       if (havebond) then
 #ifndef CLUST
 #ifndef WHAM
@@ -509,8 +530,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
index c459499..a85bff4 100644 (file)
@@ -20,7 +20,7 @@
       include "COMMON.FREE"
       include "COMMON.SBRIDGE"
       real*4 csingle(3,maxres2)
-      character*64 nazwa,bprotfile_temp
+      character*256 nazwa,bprotfile_temp
       character*3 liczba
       integer i,is,ie,j,ii,jj,k,kk,l,ll,mm,if
       integer nrec,nlines,iscor,islice
index 36ef6e6..46867c5 100644 (file)
@@ -17,7 +17,7 @@
       include 'COMMON.VAR'
       include 'COMMON.GEO'
       include 'COMMON.PROT'
-      character*64 nazwa,bprotfile_temp
+      character*256 nazwa,bprotfile_temp
       real*4 rtime,rpotE,ruconst,rt_bath,rprop(maxQ)
       double precision time
       integer iret,itmp,itraj,ntraj
index 60addc7..b64de48 100644 (file)
@@ -46,7 +46,7 @@ c      double precision tole /1.0d-1/
       double precision tt
       integer snk_p(MaxR,MaxT_h,Max_parm)
       logical lerr
-      character*64 bprotfile_temp
+      character*256 bprotfile_temp
       call opentmp(islice,ientout,bprotfile_temp)
       iii=0
       ii=0
@@ -368,7 +368,7 @@ c------------------------------------------------------------------------------
       include "COMMON.CONTACTS1"
       character*64 nazwa
       character*80 bxname,cxname
-      character*64 bprotfile_temp
+      character*256 bprotfile_temp
       character*3 liczba,licz
       character*2 licz2
       integer i,itj,ii,iii,j,k,l
index bca3771..c733c37 100644 (file)
@@ -412,7 +412,7 @@ c-------------------------------------------------------------------------------
       include "COMMON.PROTFILES"
       include "COMMON.PROT"
       include "COMMON.FREE"
-      character*64 bprotfile_temp
+      character*256 bprotfile_temp
       character*3 liczba,liczba2
       character*2 liczba1
       integer iunit,islice
@@ -472,7 +472,7 @@ c-------------------------------------------------------------------------------
       include "COMMON.SBRIDGE"
       include "COMMON.OBCINKA"
       real*4 csingle(3,maxres2)
-      character*64 nazwa,bprotfile_temp
+      character*256 nazwa,bprotfile_temp
       character*3 liczba
       character*2 liczba1
       integer i,j,ii,jj(maxslice),k,kk(maxslice),l,
index ac35de1..b397a6f 100644 (file)
@@ -23,7 +23,7 @@
       include "COMMON.SBRIDGE"
       include "COMMON.OBCINKA"
       real*4 csingle(3,maxres2)
-      character*64 nazwa,bprotfile_temp
+      character*256 nazwa,bprotfile_temp
       integer i,j,k,l,ii,jj(maxslice),kk(maxslice),ll(maxslice),
      &  mm(maxslice)
       integer iscor,islice,islice1,slice