From 7077a079d2932fdbb6cb5c28bae95bad9ff9105c Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Mon, 5 Feb 2018 13:08:33 +0100 Subject: [PATCH] corrrection in energies (no trash) --- source/cluster/CMakeLists.txt | 2 +- source/cluster/io_clust.f90 | 162 +++++++++++++++++++++++++++++++--------- source/unres/io.f90 | 34 ++++----- source/wham/conform_compar.f90 | 7 +- source/wham/enecalc.f90 | 26 ++++++- 5 files changed, 173 insertions(+), 58 deletions(-) diff --git a/source/cluster/CMakeLists.txt b/source/cluster/CMakeLists.txt index 1f39d0f..9f80cd2 100644 --- a/source/cluster/CMakeLists.txt +++ b/source/cluster/CMakeLists.txt @@ -47,7 +47,7 @@ set(UNRES_CLUSTER_WHAM_SRC0 if (Fortran_COMPILER_NAME STREQUAL "ifort") set (CMAKE_Fortran_FLAGS_RELEASE " ") set (CMAKE_Fortran_FLAGS_DEBUG "-O0 -g ") - set(FFLAGS0 "-fpp -mcmodel=medium -shared-intel -ip " ) + set(FFLAGS0 "-CB -g -fpp -mcmodel=medium -shared-intel -ip " ) elseif (Fortran_COMPILER_NAME STREQUAL "gfortran") set(FFLAGS0 "-std=legacy -mcmodel=medium " ) elseif (Fortran_COMPILER_NAME STREQUAL "pgf90") diff --git a/source/cluster/io_clust.f90 b/source/cluster/io_clust.f90 index db74bf6..7570389 100644 --- a/source/cluster/io_clust.f90 +++ b/source/cluster/io_clust.f90 @@ -4,7 +4,7 @@ use io_units ! use names use io_base !, only: ilen - use geometry_data, only: nres,c + use geometry_data, only: nres,c,nres_molec,boxxsize,boxysize,boxzsize use energy_data, only: nnt,nct,nss use control_data, only: lside implicit none @@ -865,7 +865,7 @@ subroutine add_new_cconf(jjj,jj,jj_old,icount,Next) use geometry_data, only: vbld,rad2deg,theta,phi,alph,omeg,deg2rad - use energy_data, only: itel,itype,dsc,max_ene + use energy_data, only: itel,itype,dsc,max_ene,molnum use control_data, only: symetr use geometry, only: int_from_cart1 ! implicit none @@ -882,13 +882,18 @@ ! include "COMMON.SBRIDGE" ! include "COMMON.GEO" integer :: i,j,jj,jjj,jj_old,icount,k,kk,l,ii,ib,& - nn,nn1,inan,Next,itj,chalen + nn,nn1,inan,Next,itj,chalen,mnum real(kind=8) :: etot,energia(0:max_ene) jjj=jjj+1 chalen=int((nct-nnt+2)/symetr) call int_from_cart1(.false.) do j=nnt+1,nct - if (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0) then + mnum=molnum(j) + if (mnum.eq.5) cycle + if (itype(j,mnum).eq.ntyp1_molec(mnum)) cycle + if (itype(j+1,molnum(j+1)).eq.ntyp1_molec(molnum(j+1))) cycle + + if ((vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0).and.(mnum.ne.5)) then if (j.gt.2) then if (itel(j).ne.0 .and. itel(j-1).ne.0) then write (iout,*) "Conformation",jjj,jj+1 @@ -912,8 +917,11 @@ endif enddo do j=nnt,nct - itj=itype(j) - if (itype(j).ne.10 .and. (vbld(nres+j)-dsc(iabs(itj))) & + mnum=molnum(j) + itj=itype(j,mnum) + if (mnum.eq.5) cycle + if (itype(j,1).ne.10 .and. (itype(j,mnum).ne.ntyp1_molec(mnum)) & + .and. (vbld(nres+j)-dsc(iabs(itj))) & .gt.2.0d0) then write (iout,*) "Conformation",jjj,jj+1 write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j) @@ -1248,7 +1256,8 @@ use energy_data, only: rescale_mode,distchainmax,ipot !,temp0 use control_data, only: titel,outpdb,outmol2,refstr,pdbref,& iscode,symetr,punch_dist,print_dist,nstart,nend,& - caonly,iopt,efree,lprint_cart,lprint_int + caonly,iopt,efree,lprint_cart,lprint_int,rlamb_ele,& + r_cut_ele ! implicit none ! include 'DIMENSIONS' ! include 'sizesclu.dat' @@ -1271,7 +1280,13 @@ read (INP,'(a80)') titel call card_concat(controlcard,.true.) - call readi(controlcard,'NRES',nres,0) + call readi(controlcard,'NRES',nres_molec(1),0) + call readi(controlcard,"NRES_NUCL",nres_molec(2),0) + call readi(controlcard,"NRES_CAT",nres_molec(5),0) + nres=0 + do i=1,5 + nres=nres_molec(i)+nres + enddo ! call alloc_clust_arrays allocate(rcutoff(max_cut+1)) !(max_cut+1) @@ -1296,6 +1311,13 @@ call readi(controlcard,'SYM',symetr,1) write (iout,*) 'sym', symetr call readi(controlcard,'NSTART',nstart,0) + call reada(controlcard,'BOXX',boxxsize,100.0d0) + call reada(controlcard,'BOXY',boxysize,100.0d0) + call reada(controlcard,'BOXZ',boxzsize,100.0d0) +! ions=index(controlcard,"IONS").gt.0 + call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0) + call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0) + write(iout,*) "R_CUT_ELE=",r_cut_ele call readi(controlcard,'NEND',nend,0) call reada(controlcard,'ECUT',ecut,10.0d0) call reada(controlcard,'PROB',prob_limit,0.99d0) @@ -1361,7 +1383,7 @@ real(kind=8) :: x(6*nres) !(maxvar) integer :: itype_pdb(nres) !(maxres) ! logical seq_comp - integer :: i,j,kkk + integer :: i,j,kkk,mnum ! ! Body ! @@ -1372,7 +1394,8 @@ !----------------------------- allocate(c(3,2*nres+2)) !(3,maxres2+2) maxres2=2*maxres allocate(dc(3,0:2*nres+2)) !(3,0:maxres2) - allocate(itype(nres+2)) !(maxres) + allocate(itype(nres+2,5)) !(maxres) + allocate(molnum(nres+1)) allocate(itel(nres+2)) do i=1,2*nres+2 @@ -1381,10 +1404,12 @@ dc(j,i)=0 enddo enddo + do j=1,5 do i=1,nres+2 - itype(i)=0 + itype(i,j)=0 itel(i)=0 enddo + enddo !-------------------------- ! Read weights of the subsequent energy terms. call card_concat(weightcard,.true.) @@ -1411,6 +1436,8 @@ call reada(weightcard,'SCALSCP',scalscp,1.0d0) call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0) call reada(weightcard,'DELT_CORR',delt_corr,0.5d0) + call reada(weightcard,'WCATCAT',wcatcat,0.0d0) + call reada(weightcard,'WCATPROT',wcatprot,0.0d0) call reada(weightcard,'TEMP0',temp0,300.0d0) !!! el if (index(weightcard,'SOFT').gt.0) ipot=6 ! 12/1/95 Added weight for the multi-body term WCORR @@ -1491,23 +1518,36 @@ read (inp,'(20(1x,a3))') (sequence(i),i=1,nres) endif ! Convert sequence to numeric code - do i=1,nres - itype(i)=rescode(i,sequence(i),iscode) + do i=1,nres_molec(1) + mnum=1 + molnum(i)=1 + itype(i,mnum)=rescode(i,sequence(i),iscode,mnum) + enddo + do i=nres_molec(1)+1,nres_molec(1)+nres_molec(2) + mnum=2 + molnum(i)=2 + itype(i,mnum)=rescode(i,sequence(i),iscode,mnum) + enddo + do i=nres_molec(1)+nres_molec(2)+1,nres_molec(1)+nres_molec(2)+nres_molec(5) + mnum=5 + molnum(i)=5 + itype(i,mnum)=rescode(i,sequence(i),iscode,mnum) enddo print *,nres - print '(20i4)',(itype(i),i=1,nres) + print '(20i4)',(itype(i,molnum(i)),i=1,nres) - do i=1,nres + do i=1,nres-1 #ifdef PROCOR - if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then + if (itype(i,molnum(i)).eq.ntyp1_molec(molnum(i))& + .or. itype(i+1,molnum(i+1)).eq.ntyp1_molec(molnum(i+1))) then #else - if (itype(i).eq.ntyp1) then + if (itype(i,molnum(i)).eq.ntyp1_molec(molnum(i))) then #endif itel(i)=0 #ifdef PROCOR - else if (iabs(itype(i+1)).ne.20) then + else if (iabs(itype(i+1,1)).ne.20) then #else - else if (iabs(itype(i)).ne.20) then + else if (iabs(itype(i,1)).ne.20) then #endif itel(i)=1 else @@ -1516,7 +1556,7 @@ enddo write (iout,*) "ITEL" do i=1,nres-1 - write (iout,*) i,itype(i),itel(i) + write (iout,*) i,itype(i,molnum(i)),itel(i) enddo print *,'Call Read_Bridge.' @@ -1524,8 +1564,8 @@ nnt=1 nct=nres print *,'NNT=',NNT,' NCT=',NCT - if (itype(1).eq.ntyp1) nnt=2 - if (itype(nres).eq.ntyp1) nct=nct-1 + if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2 + if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1 if (nstart.lt.nnt) nstart=nnt if (nend.gt.nct .or. nend.eq.0) nend=nct write (iout,*) "nstart",nstart," nend",nend @@ -1681,6 +1721,33 @@ write(iout,*)"po setup var" open (isidep1,file=sidepname,status="old") call getenv('SCCORPAR',sccorname) open (isccor,file=sccorname,status="old") + call getenv('THETPAR_NUCL',thetname_nucl) + open (ithep_nucl,file=thetname_nucl,status='old') + call getenv('ROTPAR_NUCL',rotname_nucl) + open (irotam_nucl,file=rotname_nucl,status='old') + call getenv('TORPAR_NUCL',torname_nucl) + open (itorp_nucl,file=torname_nucl,status='old') + call getenv('TORDPAR_NUCL',tordname_nucl) + open (itordp_nucl,file=tordname_nucl,status='old') + call getenv('SIDEPAR_NUCL',sidename_nucl) + open (isidep_nucl,file=sidename_nucl,status='old') + call getenv('SIDEPAR_SCBASE',sidename_scbase) + open (isidep_scbase,file=sidename_scbase,status='old') + call getenv('PEPPAR_PEPBASE',pepname_pepbase) + open (isidep_pepbase,file=pepname_pepbase,status='old') + call getenv('SCPAR_PHOSPH',pepname_scpho) + open (isidep_scpho,file=pepname_scpho,status='old') + call getenv('PEPPAR_PHOSPH',pepname_peppho) + open (isidep_peppho,file=pepname_peppho,status='old') + + + call getenv('LIPTRANPAR',liptranname) + open (iliptranpar,file=liptranname,status='old') + call getenv('TUBEPAR',tubename) + open (itube,file=tubename,status='old') + call getenv('IONPAR',ionname) + open (iion,file=ionname,status='old') + #ifndef OLDSCP ! ! 8/9/01 In the newest version SCp interaction constants are read from a file @@ -1696,7 +1763,7 @@ write(iout,*)"po setup var" !----------------------------------------------------------------------------- subroutine pdboutC(etot,rmsd,tytul) - use energy_data, only: ihpb,jhpb,itype + use energy_data, only: ihpb,jhpb,itype,molnum ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CONTROL' @@ -1712,16 +1779,28 @@ write(iout,*)"po setup var" 'G','H','I','J'/) integer :: ica(nres) real(kind=8) :: etot,rmsd - integer :: iatom,ichain,ires,i,j,iti - + integer :: iatom,ichain,ires,i,j,iti,mnum,mnum1,boxxxshift(3) + real(kind=8) :: boxxxx(3) write (ipdb,'(3a,1pe15.5,a,0pf7.2)') 'REMARK ',tytul(:20),& ' ENERGY ',etot,' RMS ',rmsd iatom=0 ichain=1 ires=0 + boxxxshift(1)=int(c(1,nnt)/boxxsize) +! if (c(1,nnt).lt.0) boxxxshift(1)=boxxxshift(1)-1 + boxxxshift(2)=int(c(2,nnt)/boxzsize) +! if (c(2,nnt).lt.0) boxxxshift(2)=boxxxshift(2)-1 + boxxxshift(3)=int(c(3,nnt)/boxzsize) +! if (c(3,nnt).lt.0) boxxxshift(3)=boxxxshift(3)-1 + + boxxxx(1)=boxxsize + boxxxx(2)=boxysize + boxxxx(3)=boxzsize + do i=nnt,nct - iti=itype(i) - if (iti.eq.ntyp1) then + mnum=molnum(i) + iti=itype(i,mnum) + if (iti.eq.ntyp1_molec(mnum)) then ichain=ichain+1 ires=0 write (ipdb,'(a)') 'TER' @@ -1729,27 +1808,40 @@ write(iout,*)"po setup var" ires=ires+1 iatom=iatom+1 ica(i)=iatom - write (ipdb,10) iatom,restyp(iti),chainid(ichain),& + if (mnum.eq.5) then + do j=1,3 + if ((c(j,i).gt.0).and.(c(j,nnt).lt.0)) then + c(j,i)=dmod(c(j,i),boxxxx(j))+(boxxxshift(j)-1)*boxxxx(j) + else if ((c(j,i).lt.0).and.(c(j,nnt).gt.0)) then + c(j,i)=dmod(c(j,i),boxxxx(j))+(boxxxshift(j)+1)*boxxxx(j) + else + c(j,i)=dmod(c(j,i),boxxxx(j))+(boxxxshift(j))*boxxxx(j) + endif + enddo + endif + write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),& ires,(c(j,i),j=1,3),1.0d0,tempfac(1,i) - if (iti.ne.10) then + if ((iti.ne.10).and.(mnum.ne.5)) then iatom=iatom+1 - write (ipdb,20) iatom,restyp(iti),chainid(ichain),& + write (ipdb,20) iatom,restyp(iti,mnum),chainid(ichain),& ires,(c(j,nres+i),j=1,3),1.0d0,tempfac(2,i) endif endif enddo write (ipdb,'(a)') 'TER' do i=nnt,nct-1 - if (itype(i).eq.ntyp1) cycle - if (itype(i).eq.10 .and. itype(i+1).ne.ntyp1) then + mnum=molnum(i) + mnum1=molnum(i+1) + if ((itype(i,mnum).eq.ntyp1).or.(mnum.eq.5)) cycle + if (itype(i,mnum).eq.10 .and. itype(i+1,mnum1).ne.ntyp1_molec(mnum1)) then write (ipdb,30) ica(i),ica(i+1) - else if (itype(i).ne.10 .and. itype(i+1).ne.ntyp1) then + else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum1).ne.ntyp1_molec(mnum1)) then write (ipdb,30) ica(i),ica(i+1),ica(i)+1 - else if (itype(i).ne.10 .and. itype(i+1).eq.ntyp1) then + else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) then write (ipdb,30) ica(i),ica(i)+1 endif enddo - if (itype(nct).ne.10) then + if ((itype(nct,molnum(nct)).ne.10).and.(molnum(nct).ne.5)) then write (ipdb,30) ica(nct),ica(nct)+1 endif do i=1,nss diff --git a/source/unres/io.f90 b/source/unres/io.f90 index 05f3585..ab5dc2a 100644 --- a/source/unres/io.f90 +++ b/source/unres/io.f90 @@ -769,19 +769,19 @@ call reada(weightcard,'WTURN6',wturn6,1.0D0) call reada(weightcard,'WSCCOR',wsccor,1.0D0) call reada(weightcard,'WSTRAIN',wstrain,1.0D0) - call reada(weightcard,'WVDWPP_NUCL',wvdwpp_nucl,1.0D0) - call reada(weightcard,'WELPP',welpp,1.0d0) - call reada(weightcard,'WVDWPSB',wvdwpsb,1.0d0) - call reada(weightcard,'WELPSB',welpsb,1.0D0) - call reada(weightcard,'WVDWSB',wvdwsb,1.0d0) - call reada(weightcard,'WELSB',welsb,1.0D0) - call reada(weightcard,'WBOND_NUCL',wbond_nucl,1.0D0) - call reada(weightcard,'WANG_NUCL',wang_nucl,1.0D0) - call reada(weightcard,'WSBLOC',wsbloc,1.0D0) - call reada(weightcard,'WTOR_NUCL',wtor_nucl,1.0D0) - call reada(weightcard,'WTORD_NUCL',wtor_d_nucl,1.0D0) - call reada(weightcard,'WCORR_NUCL',wcorr_nucl,1.0D0) - call reada(weightcard,'WCORR3_NUCL',wcorr3_nucl,1.0D0) + call reada(weightcard,'WVDWPP_NUCL',wvdwpp_nucl,0.0D0) + call reada(weightcard,'WELPP',welpp,0.0d0) + call reada(weightcard,'WVDWPSB',wvdwpsb,0.0d0) + call reada(weightcard,'WELPSB',welpsb,0.0D0) + call reada(weightcard,'WVDWSB',wvdwsb,0.0d0) + call reada(weightcard,'WELSB',welsb,0.0D0) + call reada(weightcard,'WBOND_NUCL',wbond_nucl,0.0D0) + call reada(weightcard,'WANG_NUCL',wang_nucl,0.0D0) + call reada(weightcard,'WSBLOC',wsbloc,0.0D0) + call reada(weightcard,'WTOR_NUCL',wtor_nucl,0.0D0) + call reada(weightcard,'WTORD_NUCL',wtor_d_nucl,0.0D0) + call reada(weightcard,'WCORR_NUCL',wcorr_nucl,0.0D0) + call reada(weightcard,'WCORR3_NUCL',wcorr3_nucl,0.0D0) call reada(weightcard,'WBOND',wbond,1.0D0) call reada(weightcard,'WTOR',wtor,1.0D0) call reada(weightcard,'WTORD',wtor_d,1.0D0) @@ -796,14 +796,14 @@ call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0) call reada(weightcard,'DELT_CORR',delt_corr,0.5d0) call reada(weightcard,'TEMP0',temp0,300.0d0) - call reada(weightcard,'WSCBASE',wscbase,1.0D0) + call reada(weightcard,'WSCBASE',wscbase,0.0D0) if (index(weightcard,'SOFT').gt.0) ipot=6 - call reada(weightcard,'WBOND_NUCL',wbond_nucl,1.0D0) + call reada(weightcard,'WBOND_NUCL',wbond_nucl,0.0D0) call reada(weightcard,'WCATCAT',wcatcat,0.0d0) call reada(weightcard,'WCATPROT',wcatprot,0.0d0) call reada(weightcard,'WPEPBASE',wpepbase,1.0d0) - call reada(weightcard,'WSCPHO',wscpho,1.0d0) - call reada(weightcard,'WPEPPHO',wpeppho,1.0d0) + call reada(weightcard,'WSCPHO',wscpho,0.0d0) + call reada(weightcard,'WPEPPHO',wpeppho,0.0d0) ! 12/1/95 Added weight for the multi-body term WCORR call reada(weightcard,'WCORRH',wcorr,1.0D0) diff --git a/source/wham/conform_compar.f90 b/source/wham/conform_compar.f90 index 0470c67..8476646 100644 --- a/source/wham/conform_compar.f90 +++ b/source/wham/conform_compar.f90 @@ -1012,7 +1012,7 @@ !---------------------------------------------------------------------------- subroutine contact(lprint,ncont,icont) - use energy_data, only: nnt,nct,itype,ipot,maxcont,sigma,sigmaii + use energy_data, only: nnt,nct,itype,ipot,maxcont,sigma,sigmaii,molnum ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' ! include 'COMMON.CHAIN' @@ -1024,6 +1024,7 @@ logical :: lprint integer :: kkk,i,j,i1,i2,it1,it2,iti,itj real(kind=8) :: rcomp + integer :: mnum,mnum2 ncont=0 kkk=3 ! print *,'nnt=',nnt,' nct=',nct @@ -3178,8 +3179,8 @@ ! include 'COMMON.CHAIN' ! include 'COMMON.NAMES' ! include 'COMMON.INTERACT' - integer :: ncont,icont(2,maxcont),isec(nres,4),nsec(nres),& - isecstr(nres) + integer :: ncont,icont(2,maxcont),isec(nres,4),nsec(nres+1),& + isecstr(nres+1) logical :: lprint,lprint_sec,not_done !el,freeres integer :: i,j,ii1,jj1,i1,j1,ij,k,ien,ist integer :: nstrand,nbeta,nhelix,iii1,jjj1,mnum diff --git a/source/wham/enecalc.f90 b/source/wham/enecalc.f90 index 597b7cc..b8fd48f 100644 --- a/source/wham/enecalc.f90 +++ b/source/wham/enecalc.f90 @@ -3,7 +3,7 @@ use io_units use wham_data ! - use geometry_data, only:nres + use geometry_data, only:nres,boxxsize,boxysize,boxzsize use energy_data ! use control_data, only:maxthetyp1 use energy, only:etotal,enerprint,rescale_weights @@ -95,7 +95,7 @@ use geometry_data, only:c,phi,theta,alph,omeg,deg2rad,anatemp,& vbld,rad2deg,dc_norm,dc,vbld_inv use io_base, only:gyrate!,briefout - use geometry, only:int_from_cart1 + use geometry, only:int_from_cart1,returnbox use io_wham, only:pdboutW use io_database, only:opentmp use conform_compar, only:qwolynes,rmsnat @@ -211,6 +211,17 @@ write(iout,*)"enecalc_ i ntot",i,ntot c(l,k+nres)=csingle(l,k+nres) enddo enddo + call returnbox + do k=1,nres + do l=1,3 + csingle(l,k)=c(l,k) + enddo + enddo + do k=nnt,nct + do l=1,3 + csingle(l,k+nres)=c(l,k+nres) + enddo + enddo anatemp= 1.0d0/(beta_h(ib,ipar)*1.987D-3) q(nQ+1,iii+1)=rmsnat(iii+1) endif @@ -1169,6 +1180,7 @@ write(iout,*)"end of store_parm" use geometry_data, only:c use io_base, only:ilen use io_wham, only:pdboutW + use geometry, only:returnbox ! implicit none ! include "DIMENSIONS" ! include "DIMENSIONS.ZSCOPT" @@ -1344,6 +1356,9 @@ write(iout,*)"end of store_parm" esccor=enetb(21,i,iparm) ! edihcnstr=enetb(20,i,iparm) edihcnstr=enetb(19,i,iparm) + ecationcation=enetb(42,i,iparm) + ecation_prot=enetb(41,i,iparm) + !#ifdef SPLITELE ! etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees & ! +wvdwpp*evdw1 & @@ -1542,6 +1557,13 @@ write(iout,*)"end of store_parm" c(k,j)=csingle(k,j) enddo enddo + call returnbox + do j=1,2*nres + do k=1,3 + csingle(k,j)=c(k,j) + enddo + enddo + eini=fdimless(i) call pdboutW(iperm(i),temper,eini,enepot(i),efree,rmsdev) enddo -- 1.7.9.5