From 56e5e2711ddd35bcd233c0e2394a6137e6e1e666 Mon Sep 17 00:00:00 2001 From: Cezary Czaplewski Date: Sat, 2 May 2020 21:04:17 +0200 Subject: [PATCH] wham & cluster PBC --- source/cluster/wham/src-HCD/DIMENSIONS | 4 +- .../wham/src-HCD/Makefile-MPICH-ifort-okeanos | 4 +- source/cluster/wham/src-HCD/geomout.F | 8 +- source/cluster/wham/src-HCD/read_constr_homology.F | 2 + source/cluster/wham/src-HCD/sizesclu.dat | 2 +- source/unres/src-HCD-5D/readpdb.F | 1 + source/wham/src-HCD/COMMON.CHAIN | 7 +- source/wham/src-HCD/DIMENSIONS | 6 +- source/wham/src-HCD/Makefile_MPICH_ifort-okeanos | 4 +- source/wham/src-HCD/enecalc1.F | 4 +- source/wham/src-HCD/molread_zs.F | 10 +- source/wham/src-HCD/oligomer.F | 183 ++++++++++++++------ source/wham/src-HCD/read_constr_homology.F | 3 +- source/wham/src-HCD/readpdb.F | 7 +- 14 files changed, 166 insertions(+), 79 deletions(-) diff --git a/source/cluster/wham/src-HCD/DIMENSIONS b/source/cluster/wham/src-HCD/DIMENSIONS index 909354c..ffd3a80 100644 --- a/source/cluster/wham/src-HCD/DIMENSIONS +++ b/source/cluster/wham/src-HCD/DIMENSIONS @@ -9,8 +9,8 @@ C Max. number of processors. parameter (maxprocs=48) C Max. number of AA residues integer maxres,maxres2 - parameter (maxres=1200) -c parameter (maxres=3300) +c parameter (maxres=1200) + parameter (maxres=5000) C Appr. max. number of interaction sites parameter (maxres2=2*maxres) C Max. number of variables diff --git a/source/cluster/wham/src-HCD/Makefile-MPICH-ifort-okeanos b/source/cluster/wham/src-HCD/Makefile-MPICH-ifort-okeanos index d1e64a7..4f7f61f 100644 --- a/source/cluster/wham/src-HCD/Makefile-MPICH-ifort-okeanos +++ b/source/cluster/wham/src-HCD/Makefile-MPICH-ifort-okeanos @@ -1,7 +1,7 @@ #INSTALL_DIR = /opt/cray/mpt/7.3.2/gni/mpich-intel/15.0 FC = ftn -#OPT = -O3 -ip -mcmodel=medium -shared-intel -dynamic -OPT = -CB -g -mcmodel=medium -shared-intel -dynamic +OPT = -O3 -ip -mcmodel=medium -shared-intel -dynamic +#OPT = -CB -g -mcmodel=medium -shared-intel -dynamic FFLAGS = ${OPT} -c -I. -Iinclude_unres -I$(INSTALL_DIR)/include LIBS = -L$(INSTALL_DIR)/lib -lmpich xdrf/libxdrf.a diff --git a/source/cluster/wham/src-HCD/geomout.F b/source/cluster/wham/src-HCD/geomout.F index 4ef656f..be43686 100644 --- a/source/cluster/wham/src-HCD/geomout.F +++ b/source/cluster/wham/src-HCD/geomout.F @@ -17,12 +17,15 @@ 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 @@ -35,6 +38,7 @@ & ires,(c(j,nres+i),j=1,3),1.0d0,tempfac(2,i) endif endif + iti_prev=iti enddo write (ipdb,'(a)') 'TER' do i=nnt,nct-1 diff --git a/source/cluster/wham/src-HCD/read_constr_homology.F b/source/cluster/wham/src-HCD/read_constr_homology.F index cac0d06..b188deb 100644 --- a/source/cluster/wham/src-HCD/read_constr_homology.F +++ b/source/cluster/wham/src-HCD/read_constr_homology.F @@ -144,8 +144,10 @@ c unres_pdb=.false. if (read2sigma) then call readpdb_template(k) + close(ipdbin) else call readpdb(out_template_coord) + close(ipdbin) endif c call readpdb diff --git a/source/cluster/wham/src-HCD/sizesclu.dat b/source/cluster/wham/src-HCD/sizesclu.dat index 7d0d666..cb8572c 100644 --- a/source/cluster/wham/src-HCD/sizesclu.dat +++ b/source/cluster/wham/src-HCD/sizesclu.dat @@ -5,7 +5,7 @@ * Max. number of conformations in the data set. * integer maxconf,maxstr_proc - PARAMETER (MAXCONF=8000) + PARAMETER (MAXCONF=12000) parameter (maxstr_proc=maxconf/2) * * Max. number of "distances" between conformations. diff --git a/source/unres/src-HCD-5D/readpdb.F b/source/unres/src-HCD-5D/readpdb.F index 1190bb1..c56b1df 100644 --- a/source/unres/src-HCD-5D/readpdb.F +++ b/source/unres/src-HCD-5D/readpdb.F @@ -30,6 +30,7 @@ C geometry. lsecondary=.false. nhfrag=0 nbfrag=0 + ires=0 do read (ipdbin,'(a80)',end=10) card if (card(:5).eq.'HELIX') then diff --git a/source/wham/src-HCD/COMMON.CHAIN b/source/wham/src-HCD/COMMON.CHAIN index 7369baa..dfffc78 100644 --- a/source/wham/src-HCD/COMMON.CHAIN +++ b/source/wham/src-HCD/COMMON.CHAIN @@ -1,6 +1,6 @@ integer nres,nres0,nsup,nstart_sup,nend_sup,nstart_seq, - & ishift_pdb,chain_length,chain_border,ichanres,tabpermchain, - & nchain ,npermchain,ireschain,iz_sc + & ishift_pdb,chain_length,chain_border,chain_border1,ichanres, + & tabpermchain,nchain ,npermchain,ireschain,iz_sc 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), @@ -11,7 +11,8 @@ & rmssing,anatemp,iz_sc,nsup, & nstart_sup,nend_sup,chain_length(maxchain),npermchain, & ireschain(maxres),tabpermchain(maxchain,maxperm), - & chain_border(2,maxchain),nchain,nstart_seq,ishift_pdb + & chain_border(2,maxchain),chain_border1(2,maxchain),nchain, + & nstart_seq,ishift_pdb double precision boxxsize,boxysize,boxzsize,enecut,sscut,sss, & sssgrad, & buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick diff --git a/source/wham/src-HCD/DIMENSIONS b/source/wham/src-HCD/DIMENSIONS index 48e0adf..4bc527f 100644 --- a/source/wham/src-HCD/DIMENSIONS +++ b/source/wham/src-HCD/DIMENSIONS @@ -14,8 +14,8 @@ c parameter (max_cg_procs=maxprocs) C Max. number of AA residues integer maxres c parameter (maxres=250) - parameter (maxres=1200) -c parameter (maxres=3300) +c parameter (maxres=1200) + parameter (maxres=5000) C Appr. max. number of interaction sites integer maxres2 parameter (maxres2=2*maxres) @@ -24,7 +24,7 @@ c Max. number of chains parameter (maxchain=6) C Max number of symetries integer maxsym,maxperm - parameter (maxsym=maxchain,maxperm=720) + parameter (maxsym=maxchain,maxperm=1200) C Max. number of variables integer maxvar parameter (maxvar=4*maxres) diff --git a/source/wham/src-HCD/Makefile_MPICH_ifort-okeanos b/source/wham/src-HCD/Makefile_MPICH_ifort-okeanos index 5280f0a..e667382 100644 --- a/source/wham/src-HCD/Makefile_MPICH_ifort-okeanos +++ b/source/wham/src-HCD/Makefile_MPICH_ifort-okeanos @@ -1,9 +1,9 @@ BIN = ~/bin FC = ftn -#OPT = -mcmodel=medium -shared-intel -O3 -dynamic +OPT = -mcmodel=medium -shared-intel -O3 -dynamic #OPT = -O3 -intel-static -mcmodel=medium #OPT = -O3 -ip -w -OPT = -g -CB -mcmodel=medium -shared-intel -dynamic +#OPT = -g -CB -mcmodel=medium -shared-intel -dynamic FFLAGS = ${OPT} -c -I. -I./include_unres -I$(INSTALL_DIR)/include LIBS = -L$(INSTALL_DIR)/lib -lmpich xdrf/libxdrf.a diff --git a/source/wham/src-HCD/enecalc1.F b/source/wham/src-HCD/enecalc1.F index f037ae8..0040e37 100644 --- a/source/wham/src-HCD/enecalc1.F +++ b/source/wham/src-HCD/enecalc1.F @@ -213,8 +213,8 @@ c call pdbout(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev) 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) diff --git a/source/wham/src-HCD/molread_zs.F b/source/wham/src-HCD/molread_zs.F index d7f586d..878e4dd 100644 --- a/source/wham/src-HCD/molread_zs.F +++ b/source/wham/src-HCD/molread_zs.F @@ -73,10 +73,18 @@ C Convert sequence to numeric code nct=nres call seq2chains(nres,itype,nchain,chain_length,chain_border, & ireschain) + chain_border1(1,1)=1 + chain_border1(2,1)=chain_border(2,1)+1 + do i=2,nchain-1 + chain_border1(1,i)=chain_border(1,i)-1 + chain_border1(2,i)=chain_border(2,i)+1 + enddo + chain_border1(1,nchain)=chain_border(1,nchain)-1 + chain_border1(2,nchain)=nres write(iout,*) "nres",nres," nchain",nchain do i=1,nchain write(iout,*)"chain",i,chain_length(i),chain_border(1,i), - & chain_border(2,i) + & 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) diff --git a/source/wham/src-HCD/oligomer.F b/source/wham/src-HCD/oligomer.F index 34b7be0..c63ea47 100644 --- a/source/wham/src-HCD/oligomer.F +++ b/source/wham/src-HCD/oligomer.F @@ -4,73 +4,142 @@ include "COMMON.CHAIN" include "COMMON.INTERACT" include "COMMON.IOUNITS" - integer i,ii,ipi,ipj,ipmin,j,jmin,k,ix,iy,iz, - & ixmin,iymin,izmin,ir_start,ir_end - integer iper(maxchain),iaux - double precision dchain,dchainmin,cmchain(3,20) - cmchain=0.0d0 - do i=1,nchain - ii=0 - do j=chain_border(1,i),chain_border(2,i) - if (itype(j).eq.ntyp1) cycle - ii=ii+1 - do k=1,3 - cmchain(k,i)=cmchain(k,i)+c(k,j) + integer i,j,k,l,ichain,jchain,ii,lshift(maxchain),lmin(maxchain), + & nrestot + double precision sumodl,sumodl_min,shift_coord(3,maxchain), + & boxsize(3),cm(3,maxchain),ccm(3),ccm_prev(3) + logical previous,change + boxsize(1)=boxxsize + boxsize(2)=boxysize + boxsize(3)=boxzsize + nrestot=0 + do ichain=1,nchain +c print *,"ichain",ichain + nrestot=nrestot+chain_length(ichain) + do j=1,3 + cm(j,ichain)=0.0d0 + enddo + do i=chain_border(1,ichain),chain_border(2,ichain) +c print *,i,c(:,i) + do j=1,3 + cm(j,ichain)=cm(j,ichain)+c(j,i) enddo enddo - do k=1,3 - cmchain(k,i)=cmchain(k,i)/ii + do j=1,3 + cm(j,ichain)=cm(j,ichain)/chain_length(ichain) enddo enddo +#ifdef DEBUG + write (iout,*) "CM before shift" do i=1,nchain - iper(i)=i + write (iout,'(i5,3f10.5)') i,(cm(j,i),j=1,3) enddo - do i=1,nchain-1 - dchainmin=1.0d10 - do j=i+1,nchain - ipi=iper(i) - ipj=iper(j) - do ix=-1,1 - do iy=-1,1 - do iz=-1,1 - dchain=(cmchain(1,ipj)-cmchain(1,ipi)+ix*boxxsize)**2+ - & (cmchain(2,ipj)-cmchain(2,ipi)+iy*boxysize)**2+ - & (cmchain(3,ipj)-cmchain(3,ipi)+iz*boxzsize)**2 -c write (iout,*) "i",i," ipi",ipi," j",j," ipj",ipj," d", -c & dsqrt(dchain)," dmin",dsqrt(dchainmin)," jmin",jmin - if (dchain.lt.dchainmin) then - dchainmin=dchain - ixmin=ix - iymin=iy - izmin=iz - jmin=j - endif - enddo +#endif DEBUG + do j=1,3 + lmin=0 + sumodl_min=1.0d10 + do i=0,3**nchain-1 + ii=i + do ichain=1,nchain + lshift(ichain)=mod(ii,3)-1 + ii=ii/3 + enddo +#ifdef DEBUG + write (iout,'(10i3)') (lshift(ichain),ichain=1,nchain) +#endif + sumodl=0.0d0 + do ichain=1,nchain + do jchain=1,ichain-1 + sumodl=sumodl+dabs(cm(j,ichain)+lshift(ichain)*boxsize(j) + & -cm(j,jchain)-lshift(jchain)*boxsize(j)) enddo enddo +#ifdef DEBUG + write(iout,*) "sumodl",sumodl," sumodl_min",sumodl_min +#endif + if (sumodl.lt.sumodl_min) then + sumodl_min=sumodl + lmin=lshift + endif enddo - if (ixmin.eq.0 .and. iymin.eq.0 .and. izmin.eq.0) cycle - ipj=iper(jmin) - cmchain(1,ipj)=cmchain(1,ipj)+ixmin*boxxsize - cmchain(2,ipj)=cmchain(2,ipj)+iymin*boxysize - cmchain(3,ipj)=cmchain(3,ipj)+izmin*boxzsize - ir_start=chain_border(1,ipj) - if (ir_start.gt.1) ir_start=ir_start-1 - ir_end=chain_border(2,ipj) - if (ir_end.lt.nres) ir_end=ir_end+1 - do k=ir_start,ir_end - c(1,k)=c(1,k)+ixmin*boxxsize - c(2,k)=c(2,k)+iymin*boxysize - c(3,k)=c(3,k)+izmin*boxzsize - c(1,k+nres)=c(1,k+nres)+ixmin*boxxsize - c(2,k+nres)=c(2,k+nres)+iymin*boxysize - c(3,k+nres)=c(3,k+nres)+izmin*boxzsize + do ichain=1,nchain + cm(j,ichain)=cm(j,ichain)+lmin(ichain)*boxsize(j) + shift_coord(j,ichain)=lmin(ichain)*boxsize(j) + enddo + enddo + do ichain=1,nchain + do j=1,3 + ccm(j)=ccm(j)+cm(j,ichain)*chain_length(ichain) enddo -c write (iout,*) "jmin",jmin," ipj",ipj, -c & " ixmin",ixmin," iymin",iymin," izmin",izmin - iaux=iper(i+1) - iper(i+1)=iper(jmin) - iper(jmin)=iaux enddo + do j=1,3 + ccm(j)=ccm(j)/nrestot + enddo +#ifdef DEBUG + write (iout,'(a,3f10.5)') "ccm ",(ccm(j),j=1,3) + write (iout,'(a,3f10.5)') "ccm_prev",(ccm_prev(j),j=1,3) + write (iout,'(a,3f10.5)') "diff ",(ccm(j)-ccm_prev(j),j=1,3) + write (iout,*) "shift_coord" + do i=1,nchain + write (iout,'(i5,3f10.5)') i,(shift_coord(j,i),j=1,3) + enddo +#endif; +#ifdef PREVIOUS +c +c Shift the system by box dimensions so that its CM be closest to the +c CM of the pevious snapshot +c + if (previous) then + do j=1,3 + change=.true. + do while (change) + if (ccm(j)-ccm_prev(j).gt.boxsize(j)/2) then + do ichain=1,nchain + shift_coord(j,ichain)=shift_coord(j,ichain)-boxsize(j) + enddo + ccm(j)=ccm(j)-boxsize(j) + else if (ccm(j)-ccm_prev(j).lt.-boxsize(j)/2) then + do ichain=1,nchain + shift_coord(j,ichain)=shift_coord(j,ichain)+boxsize(j) + enddo + ccm(j)=ccm(j)+boxsize(j) + else + change=.false. + endif + enddo + enddo + do j=1,3 + ccm_prev(j)=ccm(j) + enddo + else + previous=.true. + do j=1,3 + ccm_prev(j)=ccm(j) + enddo + endif +#else +c Simply shift the CM of the whole oligomer to the origin of the +c coordinate system + do ichain=1,nchain + do j=1,3 + shift_coord(j,ichain)=shift_coord(j,ichain)-ccm(j) + enddo + enddo +#endif +#ifdef DEBUG + write (iout,*) "shift_coord" + do i=1,nchain + write (iout,'(i5,3f10.5)') i,(shift_coord(j,i),j=1,3) + enddo +#endif + do ichain=1,nchain + do i=chain_border1(1,ichain),chain_border1(2,ichain) + do j=1,3 + c(j,i)=c(j,i)+shift_coord(j,ichain) + c(j,i+nres)=c(j,i+nres)+shift_coord(j,ichain) + enddo + enddo + enddo + return end diff --git a/source/wham/src-HCD/read_constr_homology.F b/source/wham/src-HCD/read_constr_homology.F index ab9901d..dece50f 100644 --- a/source/wham/src-HCD/read_constr_homology.F +++ b/source/wham/src-HCD/read_constr_homology.F @@ -145,8 +145,10 @@ c unres_pdb=.false. if (read2sigma) then call readpdb_template(k) + close(ipdbin) else call readpdb + close(ipdbin) endif c call readpdb @@ -163,7 +165,6 @@ c call readpdb write (iout,*) "read_constr_homology: after reading pdb file" call flush(iout) #endif - c c Distance restraints c diff --git a/source/wham/src-HCD/readpdb.F b/source/wham/src-HCD/readpdb.F index 8a85b8c..83c63ab 100644 --- a/source/wham/src-HCD/readpdb.F +++ b/source/wham/src-HCD/readpdb.F @@ -27,9 +27,10 @@ C geometry. ibeg=1 ishift1=0 sccalc=.false. + ires=0 do read (ipdbin,'(a80)',end=10) card -! write (iout,'(a)') card +c write (iout,'(a)') card if (card(:5).eq.'HELIX') then nhfrag=nhfrag+1 lsecondary=.true. @@ -58,7 +59,7 @@ C geometry. iterter(ires_old)=1 ishift1=ishift1+1 ibeg=2 -! write (iout,*) "Chain ended",ires,ishift,ires_old + write (iout,*) "Chain ended",ires,ishift,ires_old if (unres_pdb) then do j=1,3 dc(j,ires)=sccor(j,iii) @@ -119,7 +120,7 @@ c write (iout,*) "BEG ires",ires ! Start a new chain ishift=-ires_old+ires-1 !!!!! ishift1=ishift1-1 !!!!! -! write (iout,*) "New chain started",ires,ishift,ishift1,"!" + write (iout,*) "New chain started",ires,ishift,ishift1,"!" ires=ires-ishift+ishift1 ires_old=ires ibeg=0 -- 1.7.9.5