From 705644e0cbb7678faefd6fe1bc436159d38ad85d Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Mon, 29 Jan 2018 20:26:16 +0100 Subject: [PATCH] changes in wham and unres --- source/unres/MREMD.f90 | 4 +- source/unres/control.F90 | 1 + source/unres/data/names.f90 | 2 +- source/unres/energy.f90 | 12 ++- source/unres/geometry.f90 | 4 + source/unres/io.f90 | 3 +- source/unres/io_base.f90 | 3 +- source/unres/io_config.f90 | 23 ++++- source/unres/unres.f90 | 14 +-- source/wham/CMakeLists.txt | 4 +- source/wham/conform_compar.f90 | 161 +++++++++++++++++-------------- source/wham/enecalc.f90 | 29 +++--- source/wham/io_database.f90 | 32 ++++--- source/wham/io_wham.f90 | 203 +++++++++++++++++++++++++++++----------- source/wham/wham_calc.f90 | 13 ++- source/wham/wham_data.f90 | 2 +- 16 files changed, 329 insertions(+), 181 deletions(-) diff --git a/source/unres/MREMD.f90 b/source/unres/MREMD.f90 index 77be24d..522bc00 100644 --- a/source/unres/MREMD.f90 +++ b/source/unres/MREMD.f90 @@ -560,7 +560,7 @@ if (max_cache_traj_use.ne.1) & print *,itime,"processor ",me," over cache ",ntwx_cache do i=1,ntwx_cache-1 - + call returnbox totT_cache(i)=totT_cache(i+1) EK_cache(i)=EK_cache(i+1) potE_cache(i)=potE_cache(i+1) @@ -607,7 +607,7 @@ ugamma_cache(i,ntwx_cache)=ugamma(i) uscdiff_cache(i,ntwx_cache)=uscdiff(i) enddo - call returnbox +! call returnbox do i=1,nres*2 do j=1,3 c_cache(j,i,ntwx_cache)=c(j,i) diff --git a/source/unres/control.F90 b/source/unres/control.F90 index db5413f..8b908c1 100644 --- a/source/unres/control.F90 +++ b/source/unres/control.F90 @@ -1922,6 +1922,7 @@ nside=0 do i=2,nres-1 mnum=molnum(i) + write(iout,*) "i",molnum(i) #ifdef WHAM_RUN if (itype(i,1).ne.10) then #else diff --git a/source/unres/data/names.f90 b/source/unres/data/names.f90 index dc8b2e1..f3a4b12 100644 --- a/source/unres/data/names.f90 +++ b/source/unres/data/names.f90 @@ -99,7 +99,7 @@ "WELSB ","WBOND ","WANG ","WSBLOC ","WTOR ",& "WTORD ","WCORR ","WCORR3 ","WNULL ","WNULL ",& "WCATPROT ","WCATCAT ","WNULL ","WNULL ","WNULL ",& - "WSCBASE ","WPEPBASE ","WSCPHO ","WPEPPH O "/) + "WSCBASE ","WPEPBASE ","WSCPHO ","WPEPPHO "/) integer :: nprint_ene = 21 integer,dimension(n_ene) :: print_order = & diff --git a/source/unres/energy.f90 b/source/unres/energy.f90 index 4678cb5..7ab625f 100644 --- a/source/unres/energy.f90 +++ b/source/unres/energy.f90 @@ -616,6 +616,7 @@ call esb(esbloc) call multibody_hb_nucl(ecorr_nucl,ecorr3_nucl,n_corr,n_corr1) else + etors_nucl=0.0d0 estr_nucl=0.0d0 ebe_nucl=0.0d0 evdwsb=0.0d0 @@ -20419,7 +20420,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !c !c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3 !c - print *,"iatel_s_nucl,iatel_e_nucl",iatel_s_nucl,iatel_e_nucl +! print *,"iatel_s_nucl,iatel_e_nucl",iatel_s_nucl,iatel_e_nucl do i=iatel_s_nucl,iatel_e_nucl if (itype(i,2).eq.ntyp1_molec(2) .or. itype(i+1,2).eq.ntyp1_molec(2)) cycle dxi=dc(1,i) @@ -20806,8 +20807,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' evdwij=evdwij*eps2rt evdwsb=evdwsb+evdwij if (lprn) then - sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) - epsi=bb(itypi,itypj)**2/aa(itypi,itypj) + sigm=dabs(aa_nucl(itypi,itypj)/bb_nucl(itypi,itypj))**(1.0D0/6.0D0) + epsi=bb_nucl(itypi,itypj)**2/aa_nucl(itypi,itypj) write (iout,'(2(a3,i3,2x),17(0pf7.3))') & restyp(itypi,2),i,restyp(itypj,2),j, & epsi,sigm,chi1,chi2,chip1,chip2, & @@ -21764,14 +21765,17 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' r012 = r06**2 k0 = 332.0*(2.0*2.0)/80.0 itmp=0 + do i=1,4 itmp=itmp+nres_molec(i) enddo + write(iout,*) "itmp",itmp do i=itmp+1,itmp+nres_molec(5)-1 xi=c(1,i) yi=c(2,i) zi=c(3,i) + xi=mod(xi,boxxsize) if (xi.lt.0) xi=xi+boxxsize yi=mod(yi,boxysize) @@ -21790,6 +21794,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' if (yj.lt.0) yj=yj+boxysize zj=dmod(zj,boxzsize) if (zj.lt.0) zj=zj+boxzsize + write(iout,*) c(1,i),xi,xj,"xy",boxxsize dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 xj_safe=xj yj_safe=yj @@ -21847,6 +21852,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' gradcatcat(k,j)=gradcatcat(k,j)+gg(k) enddo + write(iout,*) "ecatcat",i,j, ecationcation,xj,yj,zj ecationcation=ecationcation+Evan1cat+Evan2cat+Eeleccat enddo enddo diff --git a/source/unres/geometry.f90 b/source/unres/geometry.f90 index ddb705f..2287f0a 100644 --- a/source/unres/geometry.f90 +++ b/source/unres/geometry.f90 @@ -3724,6 +3724,7 @@ allareout=1 !C change suggested by Ana -end do i=1,nres-1 + mnum=molnum(i) if ((itype(i,mnum).eq.ntyp1_molec(mnum))& .and.(itype(i+1,mnum).eq.ntyp1_molec(mnum))) then chain_end=i @@ -3952,6 +3953,9 @@ c(1,i)=dmod(c(1,i),boxxsize) c(2,i)=dmod(c(2,i),boxysize) c(3,i)=dmod(c(3,i),boxzsize) + c(1,i+nres)=dmod(c(1,i+nres),boxxsize) + c(2,i+nres)=dmod(c(2,i+nres),boxysize) + c(3,i+nres)=dmod(c(3,i+nres),boxzsize) endif enddo return diff --git a/source/unres/io.f90 b/source/unres/io.f90 index a661c7a..05f3585 100644 --- a/source/unres/io.f90 +++ b/source/unres/io.f90 @@ -1339,7 +1339,8 @@ ((c(l,k+nres),l=1,3),k=nnt,nct) write (iout,*) "Exit READ_CART" write (iout,'(8f10.5)') & - ((c(l,k),l=1,3),k=1,nres),& + ((c(l,k),l=1,3),k=1,nres) + write (iout,'(8f10.5)') & ((c(l,k+nres),l=1,3),k=nnt,nct) call int_from_cart1(.true.) write (iout,*) "Finish INT_TO_CART" diff --git a/source/unres/io_base.f90 b/source/unres/io_base.f90 index 3df6071..4be9dfe 100644 --- a/source/unres/io_base.f90 +++ b/source/unres/io_base.f90 @@ -56,7 +56,8 @@ print *,"ENTER READ" ! Read bridging residues. read (inp,*) ns - print *,"ns",ns + write(iout,*) "ns",ns + call flush(iout) if (ns.gt.0) then allocate(iss(ns)) read (inp,*) (iss(i),i=1,ns) diff --git a/source/unres/io_config.f90 b/source/unres/io_config.f90 index 7a931cc..3f4303a 100644 --- a/source/unres/io_config.f90 +++ b/source/unres/io_config.f90 @@ -2785,15 +2785,27 @@ ! write (2,*) "UNRES_PDB",unres_pdb ires=0 ires_old=0 +#ifdef WHAM_RUN + do i=1,nres + do j=1,5 + itype(i,j)=0 + enddo + enddo +#endif nres=0 iii=0 lsecondary=.false. nhfrag=0 nbfrag=0 + do j=1,5 + nres_molec(j)=0 + enddo + + !----------------------------- allocate(hfrag(2,maxres/3)) !(2,maxres/3) allocate(bfrag(4,maxres/3)) !(4,maxres/3) - + if(.not. allocated(istype)) allocate(istype(maxres)) do i=1,100000 read (ipdbin,'(a80)',end=10) card ! write (iout,'(a)') card @@ -2906,7 +2918,7 @@ ! nres_molec(molecule)=nres_molec(molecule)+1 else molecule=2 - itype(ires,molecule)=rescode(ires,res(2:4),0,molecule) + itype(ires,molecule)=rescode(ires,res(2:3),0,molecule) ! nres_molec(molecule)=nres_molec(molecule)+1 read (card(19:19),'(a1)') sugar isugar=sugarcode(sugar,ires) @@ -3029,7 +3041,8 @@ molecule=5 nres_molec(molecule)=nres_molec(molecule)+1 print *,"HERE",nres_molec(molecule) - itype(ires,molecule)=rescode(ires,res(2:4),0,molecule) + res=res(2:3)//' ' + itype(ires,molecule)=rescode(ires,res,0,molecule) read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3) endif endif !atom @@ -3280,7 +3293,9 @@ ! NOW LETS ROCK! SORTING allocate(c_temporary(3,2*nres)) allocate(itype_temporary(nres,5)) - allocate(molnum(nres+1)) + if (.not.allocated(molnum)) allocate(molnum(nres+1)) + if (.not.allocated(istype)) write(iout,*) & + "SOMETHING WRONG WITH ISTYTPE" allocate(istype_temp(nres)) itype_temporary(:,:)=0 seqalingbegin=1 diff --git a/source/unres/unres.f90 b/source/unres/unres.f90 index 3d31629..3fad3c4 100644 --- a/source/unres/unres.f90 +++ b/source/unres/unres.f90 @@ -245,7 +245,7 @@ integer :: j,k call alloc_compare_arrays - if (indpdb.eq.0) then + if ((indpdb.eq.0).and.(.not.read_cart)) then call chainbuild write(iout,*) 'Warning: Calling chainbuild' endif @@ -518,7 +518,7 @@ read (intin,'(i5)',end=1100,err=1100) iconf call read_angles(intin,*11) call geom_to_var(nvar,varia) - write(iout,*) 'Warning: Calling chainbuild' + write(iout,*) 'Warning: Calling chainbuild1' call chainbuild endif write (iout,'(a,i7)') 'Conformation #',iconf @@ -564,7 +564,7 @@ read (intin,'(i5)',end=11,err=11) iconf call read_angles(intin,*11) call geom_to_var(nvar,varia) - write(iout,*) 'Warning: Calling chainbuild' + write(iout,*) 'Warning: Calling chainbuild2' call chainbuild endif write (iout,'(a,i7)') 'Conformation #',iconf @@ -600,7 +600,7 @@ ! print *,'result received from worker ',man,' sending now' call var_to_geom(nvar,varia) - write(iout,*) 'Warning: Calling chainbuild' + write(iout,*) 'Warning: Calling chainbuild3' call chainbuild call etotal(energy_) iconf=ind(2) @@ -636,7 +636,7 @@ read (intin,'(i5)',end=1101,err=1101) iconf call read_angles(intin,*11) call geom_to_var(nvar,varia) - write(iout,*) 'Warning: Calling chainbuild' + write(iout,*) 'Warning: Calling chainbuild4' call chainbuild endif n=n+1 @@ -672,7 +672,7 @@ CG_COMM,muster,ierr) call var_to_geom(nvar,varia) - write(iout,*) 'Warning: Calling chainbuild' + write(iout,*) 'Warning: Calling chainbuild5' call chainbuild call etotal(energy_) iconf=ind(2) @@ -734,7 +734,7 @@ read (intin,'(i5)',end=11,err=11) iconf call read_angles(intin,*11) call geom_to_var(nvar,varia) - write(iout,*) 'Warning: Calling chainbuild' + write(iout,*) 'Warning: Calling chainbuild5' call chainbuild endif write (iout,'(a,i7)') 'Conformation #',iconf diff --git a/source/wham/CMakeLists.txt b/source/wham/CMakeLists.txt index dd0cb28..b9499b0 100644 --- a/source/wham/CMakeLists.txt +++ b/source/wham/CMakeLists.txt @@ -50,9 +50,9 @@ set(UNRES_WHAM_SRC0 # Set compiler flags for different sourcefiles #================================================ if (Fortran_COMPILER_NAME STREQUAL "ifort") - set (CMAKE_Fortran_FLAGS_RELEASE " ") + set (CMAKE_Fortran_FLAGS_RELEASE " -CB -g") set (CMAKE_Fortran_FLAGS_DEBUG "-O0 -g ") - set(FFLAGS0 "-fpp -mcmodel=medium -shared-intel " ) + set(FFLAGS0 "-fpp -mcmodel=medium -shared-intel " ) elseif (Fortran_COMPILER_NAME STREQUAL "gfortran") set(FFLAGS0 "-fpp -std=legacy -mcmodel=medium -g ") elseif (Fortran_COMPILER_NAME STREQUAL "pgf90") diff --git a/source/wham/conform_compar.f90 b/source/wham/conform_compar.f90 index 632dd6e..0470c67 100644 --- a/source/wham/conform_compar.f90 +++ b/source/wham/conform_compar.f90 @@ -896,7 +896,7 @@ use calc_data use geometry_data, only:c,dc,dc_norm - use energy_data, only:itype,maxcont + use energy_data, only:itype,maxcont,molnum ! implicit none ! include 'DIMENSIONS' ! include 'DIMENSIONS.ZSCOPT' @@ -909,11 +909,11 @@ ! include 'COMMON.CALC' ! include 'COMMON.CONTPAR' ! include 'COMMON.LOCAL' - integer :: ist,ien,kkk,iti,itj,itypi,itypj,i1,i2,it1,it2 + integer :: ist,ien,kkk,iti,itj,itypi,itypj,i1,i2,it1,it2,mnum,mnum2 real(kind=8) :: csc !el,dist real(kind=8),dimension(maxcont) :: cscore,omt1,omt2,omt12,& ddsc,ddla,ddlb - integer :: ncont + integer :: ncont,mhum integer,dimension(2,maxcont) :: icont real(kind=8) :: u,v,a(3),b(3),dla,dlb logical :: lprint @@ -926,18 +926,20 @@ if (lprint) then do i=1,nres mnum=molnum(i) - write (iout,110) restyp(itype(i)),i,c(1,i),c(2,i),& + write (iout,110) restyp(itype(i,mnum),mnum),i,c(1,i),c(2,i),& c(3,i),dc(1,nres+i),dc(2,nres+i),dc(3,nres+i),& dc_norm(1,nres+i),dc_norm(2,nres+i),dc_norm(3,nres+i) enddo endif 110 format (a,'(',i3,')',9f8.3) do i=ist,ien-kkk - iti=iabs(itype(i)) - if (iti.le.0 .or. iti.gt.ntyp) cycle + mnum=molnum(i) + iti=iabs(itype(i,mnum)) + if (iti.le.0 .or. iti.gt.ntyp_molec(mnum)) cycle do j=i+kkk,ien - itj=iabs(itype(j)) - if (itj.le.0 .or. itj.gt.ntyp) cycle + mnum2=molnum(j) + itj=iabs(itype(j,mnum2)) + if (itj.le.0 .or. itj.gt.ntyp_molec(mnum2)) cycle itypi=iti itypj=itj xj = c(1,nres+j)-c(1,nres+i) @@ -992,12 +994,14 @@ if (lprint) then write (iout,'(a)') 'Contact map:' do i=1,ncont + mnum=molnum(i) i1=icont(1,i) i2=icont(2,i) - it1=itype(i1) - it2=itype(i2) + it1=itype(i1,molnum(i1)) + it2=itype(i2,molnum(i2)) +! print *,"CONTACT",i1,mnum,it1,it2 write (iout,'(i3,2x,a,i4,2x,a,i4,5f8.3,3f10.5)') & - i,restyp(it1),i1,restyp(it2),i2,cscore(i),& + i,restyp(it1,mnum),i1,restyp(it2,mnum),i2,cscore(i),& sc_cutoff(iabs(it1),iabs(it2)),ddsc(i),ddla(i),ddlb(i),& omt1(i),omt2(i),omt12(i) enddo @@ -1024,9 +1028,11 @@ kkk=3 ! print *,'nnt=',nnt,' nct=',nct do i=nnt+kkk,nct - iti=iabs(itype(i)) + mnum=molnum(i) + iti=iabs(itype(i,1)) do j=nnt,i-kkk - itj=iabs(itype(j)) + mnum2=molnum(j) + itj=iabs(itype(j,1)) if (ipot.ne.4) then ! rcomp=sigmaii(iti,itj)+1.0D0 rcomp=facont*sigmaii(iti,itj) @@ -1046,12 +1052,13 @@ if (lprint) then write (iout,'(a)') 'Contact map:' do i=1,ncont + mnum=molnum(i) i1=icont(1,i) i2=icont(2,i) - it1=itype(i1) - it2=itype(i2) + it1=itype(i1,1) + it2=itype(i2,1) write (iout,'(i3,2x,a,i4,2x,a,i4)') & - i,restyp(it1),i1,restyp(it2),i2 + i,restyp(it1,mnum),i1,restyp(it2,mnum),i2 enddo endif return @@ -1090,7 +1097,7 @@ subroutine pept_cont(lprint,ncont,icont) use geometry_data, only:c - use energy_data, only:maxcont,nnt,nct,itype + use energy_data, only:maxcont,nnt,nct,itype,molnum ! implicit none ! include 'DIMENSIONS' ! include 'DIMENSIONS.ZSCOPT' @@ -1100,7 +1107,7 @@ ! include 'COMMON.FFIELD' ! include 'COMMON.NAMES' integer :: ncont,icont(2,maxcont) - integer :: i,j,k,kkk,i1,i2,it1,it2 + integer :: i,j,k,kkk,i1,i2,it1,it2,mnum logical :: lprint !el real(kind=8) :: dist real(kind=8) :: rcomp=5.5d0 @@ -1125,12 +1132,13 @@ if (lprint) then write (iout,'(a)') 'PP contact map:' do i=1,ncont + mnum=molnum(i) i1=icont(1,i) i2=icont(2,i) - it1=itype(i1) - it2=itype(i2) + it1=itype(i1,mnum) + it2=itype(i2,mnum) write (iout,'(i3,2x,a,i4,2x,a,i4)') & - i,restyp(it1),i1,restyp(it2),i2 + i,restyp(it1,mnum),i1,restyp(it2,mnum),i2 enddo endif return @@ -1141,7 +1149,7 @@ subroutine contacts_between_fragments(lprint,is,ncont,icont,& ncont_interfrag,icont_interfrag) - use energy_data, only:itype,maxcont + use energy_data, only:itype,maxcont,molnum ! include 'DIMENSIONS' ! include 'DIMENSIONS.ZSCOPT' ! include 'DIMENSIONS.COMPAR' @@ -1153,7 +1161,7 @@ integer :: icont(2,maxcont),ncont_interfrag(mmaxfrag),& icont_interfrag(2,maxcont,mmaxfrag) logical :: OK1,OK2,lprint - integer :: is,ncont,i,j,ind,nc,k,ic1,ic2,l,i1,i2,it1,it2 + integer :: is,ncont,i,j,ind,nc,k,ic1,ic2,l,i1,i2,it1,it2,mnum ! Determine the contacts that occur within a fragment and between fragments. do i=1,nfrag(1) do j=1,i @@ -1210,10 +1218,11 @@ do k=1,ncont_interfrag(ind) i1=icont_interfrag(1,k,ind) i2=icont_interfrag(2,k,ind) - it1=itype(i1) - it2=itype(i2) + mnum=molnum(i1) + it1=itype(i1,mnum) + it2=itype(i2,molnum(i2)) write (iout,'(i3,2x,a,i4,2x,a,i4)') & - i,restyp(it1),i1,restyp(it2),i2 + i,restyp(it1,mnum),i1,restyp(it2,molnum(i2)),i2 enddo enddo write (iout,*) @@ -1230,10 +1239,11 @@ do k=1,ncont_interfrag(ind) i1=icont_interfrag(1,k,ind) i2=icont_interfrag(2,k,ind) - it1=itype(i1) - it2=itype(i2) + mnum=molnum(i1) + it1=itype(i1,mnum) + it2=itype(i2,molnum(i2)) write (iout,'(i3,2x,a,i4,2x,a,i4)') & - i,restyp(it1),i1,restyp(it2),i2 + i,restyp(it1,mnum),i1,restyp(it2,molnum(i2)),i2 enddo enddo enddo @@ -1351,7 +1361,7 @@ subroutine elecont(lprint,ncont,icont,ist,ien) use geometry_data, only:c - use energy_data, only:maxcont,rpp,epp,itype,itel,vblinv,vblinv2 + use energy_data, only:maxcont,rpp,epp,itype,itel,vblinv,vblinv2,molnum ! implicit none ! include 'DIMENSIONS' ! include 'DIMENSIONS.ZSCOPT' @@ -1375,7 +1385,7 @@ real(kind=8),dimension(2,2) :: ael6c,ael3c,appc,bppc real(kind=8) :: elcutoff=-0.3d0 real(kind=8) :: elecutoff_14=-0.5d0 - integer :: ncont,icont(2,maxcont) + integer :: ncont,icont(2,maxcont),mnum real(kind=8) :: econt(maxcont) ! ! Load the constants of peptide bond - peptide bond interactions. @@ -1470,10 +1480,10 @@ do i=1,ncont i1=icont(1,i) i2=icont(2,i) - it1=itype(i1) - it2=itype(i2) + it1=itype(i1,molnum(i1)) + it2=itype(i2,molnum(i1)) write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)') & - i,restyp(it1),i1,restyp(it2),i2,econt(i) + i,restyp(it1,molnum(i1)),i1,restyp(it2,molnum(i1)),i2,econt(i) enddo endif ! For given residues keep only the contacts with the greatest energy. @@ -1554,12 +1564,14 @@ write (iout,*) write (iout,*) 'Electrostatic contacts after pruning: ' do i=1,ncont + mnum=molnum(i) i1=icont(1,i) i2=icont(2,i) - it1=itype(i1) - it2=itype(i2) + mnum=molnum(i1) + it1=itype(i1,mnum) + it2=itype(i2,molnum(i2)) write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)') & - i,restyp(it1),i1,restyp(it2),i2,econt(i) + i,restyp(it1,mnum),i1,restyp(it2,molnum(i2)),i2,econt(i) enddo endif return @@ -2362,7 +2374,7 @@ use geometry_data, only:cref,nperm use control_data, only:symetr - use energy_data, only:nnt,nct,itype + use energy_data, only:nnt,nct,itype,molnum ! implicit none ! include 'DIMENSIONS' ! include 'DIMENSIONS.ZSCOPT' @@ -2408,7 +2420,7 @@ (cref(3,jl,kkk)-cref(3,il,kkk))**2) dij=dist(il,jl) qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2) - if (itype(il).ne.10 .or. itype(jl).ne.10) then + if (itype(il,molnum(il)).ne.10 .or. itype(jl,molnum(jl)).ne.10) then nl=nl+1 d0ijCM=dsqrt( & (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ & @@ -2420,7 +2432,7 @@ qq = qq+qqij+qqijCM if (lprn) then write (iout,*) "il",il," jl",jl,& - " itype",itype(il),itype(jl) + " itype",itype(il,molnum(il)),itype(jl,molnum(jl)) write (iout,*)"d0ij",d0ij," dij",dij," d0ijCM",d0ijCM,& " dijCM",dijCM," qqij",qqij," qqijCM",qqijCM endif @@ -2449,7 +2461,7 @@ (cref(3,jl,kkk)-cref(3,il,kkk))**2) dij=dist(il,jl) qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2) - if (itype(il).ne.10 .or. itype(jl).ne.10) then + if (itype(il,molnum(il)).ne.10 .or. itype(jl,molnum(jl)).ne.10) then nl=nl+1 d0ijCM=dsqrt( & (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ & @@ -2461,7 +2473,7 @@ qq = qq+qqij+qqijCM if (lprn) then write (iout,*) "i",i," j",j," il",il," jl",jl,& - " itype",itype(il),itype(jl) + " itype",itype(il,molnum(il)),itype(jl,molnum(jl)) write (iout,*)"d0ij",d0ij," dij",dij," d0ijCM",d0ijCM,& " dijCM",dijCM," qqij",qqij," qqijCM",qqijCM endif @@ -2494,7 +2506,7 @@ (cref(3,kl,kkk)-cref(3,il,kkk))**2) dij=dist(il,kl) qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2) - if (itype(il).ne.10 .or. itype(kl).ne.10) then + if (itype(il,molnum(il)).ne.10 .or. itype(kl,molnum(kl)).ne.10) then nl=nl+1 d0ijCM=dsqrt( & (cref(1,kl+nres,kkk)-cref(1,il+nres,kkk))**2+ & @@ -2507,7 +2519,8 @@ qq = qq+qqij+qqijCM if (lprn) then write (iout,*) "i",i," j",j," k",k," l",l," il",il,& - " kl",kl," itype",itype(il),itype(kl) + " kl",kl," itype",itype(il,molnum(il)), & + itype(kl,molnum(kl)) write (iout,*) " d0ij",d0ij," dij",dij," d0ijCM",& d0ijCM," dijCM",dijCM," qqij",qqij," qqijCM",qqijCM endif @@ -2724,7 +2737,7 @@ use control_data, only:symetr use geometry_data, only:nperm,cref,c - use energy_data, only:itype + use energy_data, only:itype,molnum ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'DIMENSIONS.ZSCOPT' @@ -2748,7 +2761,7 @@ do kkk=1,nperm nnsup=0 do i=1,nres - if (itype(i).ne.ntyp1) then + if (itype(i,molnum(i)).ne.ntyp1_molec(molnum(i))) then nnsup=nnsup+1 do j=1,3 cc(j,nnsup)=c(j,i) @@ -2778,7 +2791,7 @@ subroutine define_fragments use geometry_data, only:rad2deg - use energy_data, only:itype + use energy_data, only:itype,molnum use compare_data, only:nhfrag,nbfrag,bfrag,hfrag ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' @@ -2798,7 +2811,7 @@ ! include 'COMMON.INTERACT' ! include 'COMMON.NAMES' integer :: nstrand,istrand(2,nres/2) - integer :: nhairp,ihairp(2,nres/5) + integer :: nhairp,ihairp(2,nres/5),mnum character(len=16) :: strstr(4)=reshape((/'helix','hairpin',& 'strand','strand pair'/),shape(strstr)) integer :: j,i,ii,i1,i2,i3,i4,it1,it2,it3,it4 @@ -2957,65 +2970,68 @@ write (iout,*) "The following primary fragments were found:" write (iout,*) "Helices:",nhfrag do i=1,nhfrag + mnum=molnum(i) i1=ifrag(1,1,i) i2=ifrag(2,1,i) - it1=itype(i1) - it2=itype(i2) + it1=itype(i1,mnum) + it2=itype(i2,mnum) write (iout,'(i3,2x,a,i4,2x,a,i4)') & - i,restyp(it1),i1,restyp(it2),i2 + i,restyp(it1,mnum),i1,restyp(it2,mnum),i2 enddo write (iout,*) "Hairpins:",nhairp do i=nhfrag+1,nhfrag+nhairp i1=ifrag(1,1,i) i2=ifrag(2,1,i) - it1=itype(i1) - it2=itype(i2) + it1=itype(i1,molnum(i1)) + it2=itype(i2,molnum(i2)) write (iout,'(i3,2x,a,i4,2x,a,i4,2x)') & - i,restyp(it1),i1,restyp(it2),i2 + i,restyp(it1,molnum(i1)),i1,restyp(it2,molnum(i2)),i2 enddo write (iout,*) "Far strand pairs:",nbfrag do i=nhfrag+nhairp+1,nhfrag+nhairp+nbfrag i1=ifrag(1,1,i) i2=ifrag(2,1,i) - it1=itype(i1) - it2=itype(i2) + it1=itype(i1,molnum(i1)) + it2=itype(i2,molnum(i1)) i3=ifrag(1,2,i) i4=ifrag(2,2,i) - it3=itype(i3) - it4=itype(i4) + it3=itype(i3,molnum(i3)) + it4=itype(i4,molnum(i4)) write (iout,'(i3,2x,a,i4,2x,a,i4," and ",a,i4,2x,a,i4)') & - i,restyp(it1),i1,restyp(it2),i2,& - restyp(it3),i3,restyp(it4),i4 + i,restyp(it1,molnum(i1)),i1,restyp(it2,molnum(i2)),i2,& + restyp(it3,molnum(i3)),i3,restyp(it4,molnum(i4)),i4 enddo write (iout,*) "Strands:",nstrand do i=nhfrag+nhairp+nbfrag+1,nfrag(1) + mnum=molnum(i) i1=ifrag(1,1,i) i2=ifrag(2,1,i) - it1=itype(i1) - it2=itype(i2) + it1=itype(i1,mnum) + it2=itype(i2,mnum) write (iout,'(i3,2x,a,i4,2x,a,i4)') & - i,restyp(it1),i1,restyp(it2),i2 + i,restyp(it1,mnum),i1,restyp(it2,mnum),i2 enddo call imysort(nfrag(1),2,maxpiece,ifrag(1,1,1),npiece(1,1),& istruct(1),n_shift(1,1,1),ang_cut(1),ang_cut1(1),frac_min(1),& nc_fragm(1,1),nc_req_setf(1,1)) write (iout,*) "Fragments after sorting:" do i=1,nfrag(1) + mnum=molnum(i) i1=ifrag(1,1,i) i2=ifrag(2,1,i) - it1=itype(i1) - it2=itype(i2) + it1=itype(i1,mnum) + it2=itype(i2,mnum) write (iout,'(i3,2x,a,i4,2x,a,i4,$)') & - i,restyp(it1),i1,restyp(it2),i2 + i,restyp(it1,molnum(it1)),i1,restyp(it2,molnum(it2)),i2 if (npiece(i,1).eq.1) then write (iout,'(2x,a)') strstr(istruct(i)) else i1=ifrag(1,2,i) i2=ifrag(2,2,i) - it1=itype(i1) - it2=itype(i2) + it1=itype(i1,mnum) + it2=itype(i2,mnum) write (iout,'(2x,a,i4,2x,a,i4,2x,a)') & - restyp(it1),i1,restyp(it2),i2,strstr(istruct(i)) + restyp(it1,molnum(it1)),i1,restyp(it2,molnum(it2)),i2,strstr(istruct(i)) endif enddo return @@ -3149,7 +3165,7 @@ subroutine secondary2(lprint,lprint_sec,ncont,icont,isecstr) use geometry_data, only:anatemp,rad2deg,phi,nstart_sup,nend_sup - use energy_data, only:itype,maxcont + use energy_data, only:itype,maxcont,molnum use compare_data, only:bfrag,hfrag,nbfrag,nhfrag use compare, only:freeres ! implicit real*8 (a-h,o-z) @@ -3166,7 +3182,7 @@ isecstr(nres) 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 + integer :: nstrand,nbeta,nhelix,iii1,jjj1,mnum real(kind=8) :: p1,p2 !rel external freeres character(len=1) :: csec(0:2)=reshape((/'-','E','H'/),shape(csec)) @@ -3478,11 +3494,12 @@ write (iout,*) write (iout,*) "Secondary structure" do i=1,nres,80 + mnum=molnum(i) ist=i ien=min0(i+79,nres) write (iout,*) write (iout,'(8(7x,i3))') (k,k=ist+9,ien,10) - write (iout,'(80a1)') (onelet(itype(k)),k=ist,ien) + write (iout,'(80a1)') (onelet(itype(k,mnum)),k=ist,ien) write (iout,'(80a1)') (csec(isecstr(k)),k=ist,ien) enddo write (iout,*) diff --git a/source/wham/enecalc.f90 b/source/wham/enecalc.f90 index ecbce9d..2c432f5 100644 --- a/source/wham/enecalc.f90 +++ b/source/wham/enecalc.f90 @@ -318,7 +318,7 @@ write(iout,*)"enecalc_ i ntot",i,ntot write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct) ! call enerprint(energia(0),fT) call enerprint(energia(0)) - write (iout,'(2i5,21f8.2)') i,iparm,(energia(k),k=1,21) + write (iout,'(2i5,21f8.2)') i,iparm,(energia(k),k=1,49) write (iout,*) "ftors",ftors !el call briefout(i,energia(0)) temp=1.0d0/(beta_h(ib,ipar)*1.987D-3) @@ -393,7 +393,7 @@ write(iout,*)"enecalc_ i ntot",i,ntot endif endif potE(iii+1,iparm)=energia(0) - do k=1,21 + do k=1,49 enetb(k,iii+1,iparm)=energia(k) enddo ! write (iout,'(2i5,21f8.2)') "debug",k,iii+1,(enetb(k,iii+1,iparm),k=1,21) @@ -530,14 +530,17 @@ write(iout,*)"enecalc_ i ntot",i,ntot #ifdef MPI integer :: IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE) #endif - integer :: j,k,l,ii,itj,iprint + integer :: j,k,l,ii,itj,iprint,mnum,mnum1 if (.not.check_conf) then conf_check=.true. return endif call int_from_cart1(.false.) do j=nnt+1,nct - if (itype(j-1).ne.ntyp1 .and. itype(j).ne.ntyp1 .and. & + mnum=molnum(j) + mnum1=molnum(j-1) + if (mnum.ne.1) cycle + if (itype(j-1,mnum1).ne.ntyp1_molec(mnum1) .and. itype(j,mnum).ne.ntyp1_molec(mnum) .and. & (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0)) then if (iprint.gt.0) & write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),& @@ -561,8 +564,10 @@ write(iout,*)"enecalc_ i ntot",i,ntot endif enddo do j=nnt,nct - itj=itype(j) - if (itype(j).ne.10 .and.itype(j).ne.ntyp1 .and. & + mnum=molnum(j) + if (mnum.gt.3) cycle + itj=itype(j,mnum) + if (itype(j,mnum).ne.10 .and.itype(j,mnum).ne.ntyp1 .and. & (vbld(nres+j)-dsc(iabs(itj))).gt.2.0d0) then if (iprint.gt.0) & write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j),& @@ -841,8 +846,8 @@ write(iout,*)"enecalc_ i ntot",i,ntot ! Store sidechain parameters do i=1,ntyp do j=1,ntyp - aa_all(j,i,iparm)=aa(j,i) - bb_all(j,i,iparm)=bb(j,i) + aa_all(j,i,iparm)=aa_aq(j,i) + bb_all(j,i,iparm)=bb_aq(j,i) r0_all(j,i,iparm)=r0(j,i) sigma_all(j,i,iparm)=sigma(j,i) chi_all(j,i,iparm)=chi(j,i) @@ -1112,8 +1117,8 @@ write(iout,*)"end of store_parm" ! Restore sidechain parameters do i=1,ntyp do j=1,ntyp - aa(j,i)=aa_all(j,i,iparm) - bb(j,i)=bb_all(j,i,iparm) + aa_aq(j,i)=aa_all(j,i,iparm) + bb_aq(j,i)=bb_all(j,i,iparm) r0(j,i)=r0_all(j,i,iparm) sigma(j,i)=sigma_all(j,i,iparm) chi(j,i)=chi_all(j,i,iparm) @@ -1369,7 +1374,7 @@ write(iout,*)"end of store_parm" +wturn3*eello_turn3 & +wturn6*eello_turn6+wel_loc*eel_loc & +edihcnstr+wtor_d*etors_d+wsccor*esccor & - +wbond*estr + +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation #else etot=wsc*evdw+wscp*evdw2 & +welec*(ees+evdw1) & @@ -1379,7 +1384,7 @@ write(iout,*)"end of store_parm" +wturn3*eello_turn3 & +wturn6*eello_turn6+wel_loc*eel_loc+edihcnstr & +wtor_d*etors_d+wsccor*esccor & - +wbond*estr + +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation #endif #ifdef MPI diff --git a/source/wham/io_database.f90 b/source/wham/io_database.f90 index c501d9f..e5ab5bc 100644 --- a/source/wham/io_database.f90 +++ b/source/wham/io_database.f90 @@ -394,7 +394,7 @@ write(iout,*) "end of read database" !-------------------------------------------------------------------------------- subroutine cxread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm,*) -#define DEBUG +!#define DEBUG #ifdef DEBUG use geometry, only:int_from_cart1 use geometry_data, only:vbld,rad2deg,theta,phi,alph,omeg @@ -469,7 +469,7 @@ write(iout,*) "end of read database" call xdrffloat_(ixdrf, rtime, iret) print *,"rtime",rtime," iret",iret !d call xdrffloat_(ixdrf, rpotE, iret) - write (iout,*) "rpotE",rpotE," iret",iret !d +! write (iout,*) "rpotE",rpotE," iret",iret !d call flush(iout) call xdrffloat_(ixdrf, ruconst, iret) call xdrffloat_(ixdrf, rt_bath, iret) @@ -487,7 +487,7 @@ write(iout,*) "end of read database" #else call xdrffloat(ixdrf, rtime, iret) call xdrffloat(ixdrf, rpotE, iret) - write (iout,*) "rpotE",rpotE," iret",iret !d +! write (iout,*) "rpotE",rpotE," iret",iret !d call flush(iout) call xdrffloat(ixdrf, ruconst, iret) call xdrffloat(ixdrf, rt_bath, iret) @@ -497,7 +497,7 @@ write(iout,*) "end of read database" call xdrfint(ixdrf, jhpb(j), iret) enddo call xdrfint(ixdrf, nprop, iret) - write (iout,*) "nprop",nprop !d +! write (iout,*) "nprop",nprop !d if (it.gt.0 .and. nprop.ne.nprop_prev) then write (iout,*) "Warning previous nprop",nprop_prev,& " current",nprop @@ -514,7 +514,7 @@ write(iout,*) "end of read database" #endif if (iret.eq.0) exit itraj=mod(it,totraj(iR,iparm)) -#define DEBUG +!#define DEBUG #ifdef DEBUG write (iout,*) "ii",ii," itraj",itraj," it",it #endif @@ -542,7 +542,7 @@ write(iout,*) "end of read database" #ifdef DEBUG write (iout,'(10f8.3)') ((xoord(j,i),j=1,3),i=1,2*nres+2) #endif -#undef DEBUG +!#undef DEBUG if (iret.eq.0) exit if (itmp .ne. nres + nct - nnt + 1) then write (iout,*) "Error: inconsistent sizes",itmp,nres+nct-nnt+1 @@ -551,10 +551,10 @@ write(iout,*) "end of read database" endif time=rtime - write (iout,*) "calling slice" !d +! write (iout,*) "calling slice" !d call flush(iout) !d islice=slice(nstep(itraj),time,is,ie,ts,te) - write (iout,*) "islice",islice !d +! write (iout,*) "islice",islice !d call flush(iout) !d do i=1,nres @@ -704,7 +704,7 @@ write(iout,*) "end of read database" " conformations stored so far, slice",islice enddo call flush(iout) -#undef DEBUG +!#undef DEBUG return end subroutine cxread !-------------------------------------------------------------------------------- @@ -1371,7 +1371,7 @@ write(iout,*) "end of read database" use names, only:ntyp1 use geometry_data - use energy_data, only:itype,dsc + use energy_data, only:itype,dsc,molnum use geometry, only:int_from_cart1 ! use ! include "DIMENSIONS" @@ -1403,14 +1403,16 @@ write(iout,*) "end of read database" include "mpif.h" integer :: IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE) #endif - integer :: j,k,l,ii,itj,iprint + integer :: j,k,l,ii,itj,iprint,mnum if (.not. check_conf) then conf_check=.true. return endif call int_from_cart1(.false.) do j=nnt+1,nct - if (itype(j-1).ne.ntyp1 .and. itype(j).ne.ntyp1 .and. & + mnum=molnum(j) + if (mnum.eq.5) cycle + if (itype(j-1,mnum).ne.ntyp1 .and. itype(j,mnum).ne.ntyp1_molec(mnum) .and. & (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0)) then if (iprint.gt.0) & write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),& @@ -1434,8 +1436,10 @@ write(iout,*) "end of read database" endif enddo do j=nnt,nct - itj=itype(j) - if (itype(j).ne.10 .and.itype(j).ne.ntyp1 .and. & + mnum=molnum(j) + if (mnum.eq.5) cycle + itj=itype(j,mnum) + if (itype(j,mnum).ne.10 .and.itype(j,mnum).ne.ntyp1_molec(mnum) .and. & (vbld(nres+j)-dsc(iabs(itj))).gt.2.0d0) then if (iprint.gt.0) & write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j),& diff --git a/source/wham/io_wham.f90 b/source/wham/io_wham.f90 index 399dc1a..eab5a49 100644 --- a/source/wham/io_wham.f90 +++ b/source/wham/io_wham.f90 @@ -131,7 +131,7 @@ ! Read molecular data. ! use energy_data - use geometry_data, only:nres,deg2rad,c,dc + use geometry_data, only:nres,deg2rad,c,dc,nres_molec use control_data, only:iscode use control, only:rescode,setup_var,init_int_table use geometry, only:alloc_geo_arrays @@ -154,10 +154,17 @@ !el integer :: rescode !el real(kind=8) :: x(maxvar) character(len=320) :: controlcard !,ucase - integer,dimension(nres) :: itype_pdb !(maxres) - integer :: i,j,i1,i2,it1,it2 + integer,dimension(nres,5) :: itype_pdb !(maxres) + integer :: i,j,i1,i2,it1,it2,mnum real(kind=8) :: scalscp !el logical :: seq_comp + write(iout,*) "START MOLREAD" + call flush(iout) + do i=1,5 + nres_molec(i)=0 + print *,"nres_molec, initial",nres_molec(i) + enddo + call card_concat(controlcard,.true.) call reada(controlcard,'SCAL14',scal14,0.4d0) call reada(controlcard,'SCALSCP',scalscp,1.0d0) @@ -165,8 +172,14 @@ call reada(controlcard,'TEMP0',temp0,300.0d0) !el call reada(controlcard,'DELT_CORR',delt_corr,0.5d0) r0_corr=cutoff_corr-delt_corr - call readi(controlcard,"NRES",nres,0) - allocate(sequence(nres+1)) + 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 +! allocate(sequence(nres+1)) !el znamy juz ilosc reszt wiec mozna zaalokowac tablice do liczenia enerii call alloc_geo_arrays call alloc_ener_arrays @@ -174,8 +187,9 @@ !---------------------------- allocate(c(3,2*nres+2)) allocate(dc(3,0:2*nres+2)) - allocate(itype(nres+2)) + allocate(itype(nres+2,5)) allocate(itel(nres+2)) + if (.not. allocated(molnum)) allocate(molnum(nres+2)) ! ! Zero out tableis. do i=1,2*nres+2 @@ -185,11 +199,15 @@ enddo enddo do i=1,nres+2 - itype(i)=0 + molnum(i)=0 + do j=1,5 + itype(i,j)=0 + enddo itel(i)=0 enddo !-------------------------- ! + allocate(sequence(nres+1)) iscode=index(controlcard,"ONE_LETTER") if (nres.le.0) then write (iout,*) "Error: no residues in molecule" @@ -199,6 +217,7 @@ write (iout,*) "Error: too many residues",nres,maxres endif write(iout,*) 'nres=',nres +! HERE F**** CHANGE ! Read sequence of the protein if (iscode.gt.0) then read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres) @@ -206,22 +225,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 + write (iout,*) "Numeric code:" - write (iout,'(20i4)') (itype(i),i=1,nres) + write (iout,'(20i4)') (itype(i,molnum(i)),i=1,nres) do i=1,nres-1 + mnum=molnum(i) #ifdef PROCOR - if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then + if (itype(i,mnum).eq.ntyp1_molec(mnum) .or. itype(i+1,mnum).eq.ntyp1_molec(mnum)) then #else - if (itype(i).eq.ntyp1) then + if (itype(i,mnum).eq.ntyp1_molec(mnum)) then #endif itel(i)=0 #ifdef PROCOR - else if (iabs(itype(i+1)).ne.20) then + else if (iabs(itype(i+1,mnum)).ne.20) then #else - else if (iabs(itype(i)).ne.20) then + else if (iabs(itype(i,mnum)).ne.20) then #endif itel(i)=1 else @@ -230,8 +263,10 @@ enddo write (iout,*) "ITEL" do i=1,nres-1 - write (iout,*) i,itype(i),itel(i) + mnum=molnum(i) + write (iout,*) i,itype(i,mnum),itel(i) enddo + write(iout,*) call read_bridge if (with_dihed_constr) then @@ -256,8 +291,8 @@ nnt=1 nct=nres - 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 write(iout,*) 'NNT=',NNT,' NCT=',NCT call setup_var call init_int_table @@ -267,12 +302,13 @@ write (iout,'(20i4)') (iss(i),i=1,ns) write (iout,'(/a/)') 'Pre-formed links are:' do i=1,nss + mnum=molnum(i) i1=ihpb(i)-nres i2=jhpb(i)-nres - it1=itype(i1) - it2=itype(i2) + it1=itype(i1,mnum) + it2=itype(i2,mnum) write (iout,'(2a,i3,3a,i3,a,3f10.3)') & - restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',& + restyp(it1,molnum(i1)),'(',i1,') -- ',restyp(it2,molnum(i2)),'(',i2,')',& dhpb(i),ebr,forcon(i) enddo endif @@ -339,12 +375,12 @@ character(len=16) :: key integer :: iparm !el real(kind=8) :: ip,mp - real(kind=8) :: dwa16,akl,si,rri,epsij,rrij,sigeps,sigt1sq,& + real(kind=8) :: dwa16,akl,si,rri,epsij,rrij,sigeps,sigt1sq,epsijlip,& sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm real(kind=8) :: v0ij,v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,rjunk,& res1 integer :: i,j,ichir1,ichir2,k,l,m,kk,ii,mm,junk,lll,ll,llll,n - integer :: nlobi,iblock,maxinter,iscprol + integer :: nlobi,iblock,maxinter,iscprol,ncatprotparm ! ! Body ! @@ -391,6 +427,8 @@ allocate(ww(max_eneW)) wvdwpp=ww(16) wbond=ww(18) wsccor=ww(19) + wcatcat=ww(44) + wcatprot=ww(43) endif ! @@ -415,6 +453,9 @@ allocate(ww(max_eneW)) weights(17)=wbond weights(18)=0 !scal14 ! weights(21)=wsccor + weights(43)=wcatprot + weights(44)=wcatcat + ! el-------- call card_concat(controlcard,.false.) @@ -452,6 +493,9 @@ allocate(ww(max_eneW)) call reads(controlcard,"SCPPAR",scpname_t,scpname) open (iscpp,file=scpname_t,status='old') rewind(iscpp) + call getenv_loc('IONPAR',ionname) + open (iion,file=ionname,status='old') + rewind(iion) write (iout,*) "Parameter set:",iparm write (iout,*) "Energy-term weights:" do i=1,n_eneW @@ -507,7 +551,7 @@ allocate(ww(max_eneW)) endif enddo #else - read (ibond,*) ijunk,vbldp0,akp,rjunk + read (ibond,*) ijunk,vbldp0,vbldpDUM,akp,rjunk do i=1,ntyp read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),& j=1,nbondterm(i)) @@ -525,7 +569,7 @@ allocate(ww(max_eneW)) 'inertia','Pstok' write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0 do i=1,ntyp - write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),& + write (iout,'(a10,i3,6f10.5)') restyp(i,molnum(i)),nbondterm(i),& vbldsc0(1,i),aksc(1,i),abond0(1,i) do j=2,nbondterm(i) write (iout,'(13x,3f10.5)') & @@ -533,6 +577,23 @@ allocate(ww(max_eneW)) enddo enddo endif + if (.not. allocated(msc)) allocate(msc(ntyp1,5)) + if (.not. allocated(restok)) allocate(restok(ntyp1,5)) + + do i=1,ntyp_molec(5) + read(iion,*) msc(i,5),restok(i,5) + print *,msc(i,5),restok(i,5) + enddo + ip(5)=0.2 +! isc(5)=0.2 + read (iion,*) ncatprotparm + allocate(catprm(ncatprotparm,4)) + do k=1,4 + read (iion,*) (catprm(i,k),i=1,ncatprotparm) + enddo + print *, catprm + + !---------------------------------------------------- allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp)) allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp) @@ -1515,6 +1576,7 @@ enddo allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2) allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp) allocate(chip(ntyp1),alp(ntyp1)) !(ntyp) + allocate(epslip(ntyp,ntyp)) do i=1,ntyp do j=1,ntyp augm(i,j)=0.0D0 @@ -1548,7 +1610,7 @@ enddo call printmat(ntyp,ntyp,ntyp,iout,restyp,eps) write (iout,'(/a)') 'One-body parameters:' write (iout,'(a,4x,a)') 'residue','sigma' - write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp) + write (iout,'(a3,6x,f10.5)') (restyp(i,molnum(i)),sigma0(i),i=1,ntyp) endif ! goto 50 !----------------------- LJK potential -------------------------------- @@ -1562,7 +1624,7 @@ enddo call printmat(ntyp,ntyp,ntyp,iout,restyp,eps) write (iout,'(/a)') 'One-body parameters:' write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 ' - write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),& + write (iout,'(a3,6x,2f10.5)') (restyp(i,molnum(i)),sigma0(i),rr0(i),& i=1,ntyp) endif ! goto 50 @@ -1576,6 +1638,9 @@ enddo read (isidep,*)(sigii(i),i=1,ntyp) read (isidep,*)(chip(i),i=1,ntyp) read (isidep,*)(alp(i),i=1,ntyp) + do i=1,ntyp + read (isidep,*)(epslip(i,j),j=i,ntyp) + enddo ! For the GB potential convert sigma'**2 into chi' if (ipot.eq.4) then do i=1,ntyp @@ -1589,7 +1654,7 @@ enddo write (iout,'(/a)') 'One-body parameters:' write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',& ' chip ',' alph ' - write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),& + write (iout,'(a3,6x,4f10.5)') (restyp(i,molnum(i)),sigma0(i),sigii(i),& chip(i),alp(i),i=1,ntyp) endif ! goto 50 @@ -1606,7 +1671,7 @@ enddo write (iout,'(/a)') 'One-body parameters:' write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',& 's||/s_|_^2',' chip ',' alph ' - write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),& + write (iout,'(a3,6x,5f10.5)') (restyp(i,molnum(i)),sigma0(i),rr0(i),& sigii(i),chip(i),alp(i),i=1,ntyp) endif case default @@ -1621,12 +1686,16 @@ enddo ! Calculate the "working" parameters of SC interactions. !el from module energy - COMMON.INTERACT------- - allocate(aa(ntyp1,ntyp1),bb(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp) +! allocate(aa(ntyp1,ntyp1),bb(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp) + allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp) + allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp) allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1) do i=1,ntyp1 do j=1,ntyp1 - aa(i,j)=0.0D0 - bb(i,j)=0.0D0 + aa_aq(i,j)=0.0D0 + bb_aq(i,j)=0.0D0 + aa_lip(i,j)=0.0d0 + bb_lip(i,j)=0.0d0 chi(i,j)=0.0D0 sigma(i,j)=0.0D0 r0(i,j)=0.0D0 @@ -1637,6 +1706,7 @@ enddo do i=2,ntyp do j=1,i-1 eps(i,j)=eps(j,i) + epslip(i,j)=epslip(j,i) enddo enddo do i=1,ntyp @@ -1665,10 +1735,17 @@ enddo epsij=eps(i,j) sigeps=dsign(1.0D0,epsij) epsij=dabs(epsij) - aa(i,j)=epsij*rrij*rrij - bb(i,j)=-sigeps*epsij*rrij - aa(j,i)=aa(i,j) - bb(j,i)=bb(i,j) + aa_aq(i,j)=epsij*rrij*rrij + bb_aq(i,j)=-sigeps*epsij*rrij + aa_aq(j,i)=aa_aq(i,j) + bb_aq(j,i)=bb_aq(i,j) + epsijlip=epslip(i,j) + sigeps=dsign(1.0D0,epsijlip) + epsijlip=dabs(epsijlip) + aa_lip(i,j)=epsijlip*rrij*rrij + bb_lip(i,j)=-sigeps*epsijlip*rrij + aa_lip(j,i)=aa_lip(i,j) + bb_lip(j,i)=bb_lip(i,j) if (ipot.gt.2) then sigt1sq=sigma0(i)**2 sigt2sq=sigma0(j)**2 @@ -1701,7 +1778,7 @@ enddo endif if (lprint) then write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') & - restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),& + restyp(i,molnum(i)),restyp(j,molnum(i)),aa_aq(i,j),bb_aq(i,j),augm(i,j),& sigma(i,j),r0(i,j),chi(i,j),chi(j,i) endif enddo @@ -1874,8 +1951,9 @@ enddo !----------------------------------------------------------------------------- subroutine read_general_data(*) - use control_data, only:indpdb,symetr + use control_data, only:indpdb,symetr,r_cut_ele,rlamb_ele,ions use energy_data, only:distchainmax + use geometry_data, only:boxxsize,boxysize,boxzsize ! implicit none ! include "DIMENSIONS" ! include "DIMENSIONS.ZSCOPT" @@ -1952,6 +2030,13 @@ enddo call reada(controlcard,"DELTRMS",deltrms,5.0d-2) call reada(controlcard,"DELTRGY",deltrgy,5.0d-2) call readi(controlcard,"RESCALE",rescale_modeW,1) + 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 check_conf=index(controlcard,"NO_CHECK_CONF").eq.0 call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0) call readi(controlcard,'SYM',symetr,1) @@ -2494,7 +2579,7 @@ enddo use control_data, only:pdbref use geometry_data, only:nres,cref,c,dc,nsup,dc_norm,nend_sup,& nstart_sup,nstart_seq,nperm,nres0 - use energy_data, only:nct,nnt,icont_ref,ncont_ref,itype + use energy_data, only:nct,nnt,icont_ref,ncont_ref,itype,molnum use compare, only:seq_comp !,contact,elecont use geometry, only:chainbuild,dist use io_config, only:readpdb @@ -2522,9 +2607,9 @@ enddo character(len=4) :: sequence(nres) !el integer rescode !el real(kind=8) :: x(maxvar) - integer :: itype_pdb(nres) + integer :: itype_pdb(nres,5) !el logical seq_comp - integer :: i,j,k,nres_pdb,iaux + integer :: i,j,k,nres_pdb,iaux,mnum real(kind=8) :: ddsc !el,dist integer :: kkk !,ilen !el external ilen @@ -2541,15 +2626,16 @@ enddo return 1 34 continue do i=1,nres - itype_pdb(i)=itype(i) + mnum=molnum(i) + itype_pdb(i,mnum)=itype(i,mnum) enddo call readpdb do i=1,nres - iaux=itype_pdb(i) - itype_pdb(i)=itype(i) - itype(i)=iaux + iaux=itype_pdb(i,mnum) + itype_pdb(i,mnum)=itype(i,mnum) + itype(i,mnum)=iaux enddo close (ipdbin) do kkk=1,nperm @@ -2558,7 +2644,7 @@ enddo nstart_seq=nnt if (nsup.le.(nct-nnt+1)) then do i=0,nct-nnt+1-nsup - if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),& + if (seq_comp(itype(nnt+i,molnum(nnt+i)),itype_pdb(nstart_sup,molnum(nstart_sup)),& nsup)) then do j=nnt+nsup-1,nnt,-1 do k=1,3 @@ -2590,7 +2676,7 @@ enddo return 1 else do i=0,nsup-(nct-nnt+1) - if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),& + if (seq_comp(itype(nnt,molnum(nnt)),itype_pdb(nstart_sup+i,molnum(nstart_sup+i)),& nct-nnt+1)) & then nstart_sup=nstart_sup+i @@ -2628,10 +2714,11 @@ enddo enddo enddo do i=1,nres + mnum=molnum(i) do j=1,3 dc(j,nres+i)=cref(j,nres+i,kkk)-cref(j,i,kkk) enddo - if (itype(i).ne.10) then + if (itype(i,mnum).ne.10) then ddsc = dist(i,nres+i) do j=1,3 dc_norm(j,nres+i)=dc(j,nres+i)/ddsc @@ -2671,7 +2758,7 @@ enddo subroutine pdboutW(ii,temp,efree,etot,entropy,rmsdev) use geometry_data, only:nres,c - use energy_data, only:nss,nnt,nct,ihpb,jhpb,itype + use energy_data, only:nss,nnt,nct,ihpb,jhpb,itype,molnum ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'DIMENSIONS.ZSCOPT' @@ -2686,7 +2773,7 @@ enddo 'D','E','F','G','H','I','J'/),shape(chainid)) integer,dimension(nres) :: ica !(maxres) real(kind=8) :: temp,efree,etot,entropy,rmsdev - integer :: ii,i,j,iti,ires,iatom,ichain + integer :: ii,i,j,iti,ires,iatom,ichain,mnum write(ipdb,'("REMARK CONF",i8," TEMPERATURE",f7.1," RMS",0pf7.2)')& ii,temp,rmsdev write (ipdb,'("REMARK DIMENSIONLESS FREE ENERGY",1pe15.5)') & @@ -2697,7 +2784,8 @@ enddo ichain=1 ires=0 do i=nnt,nct - iti=itype(i) + mnum=molnum(i) + iti=itype(i,mnum) if (iti.eq.ntyp1) then ichain=ichain+1 ires=0 @@ -2706,27 +2794,28 @@ enddo ires=ires+1 iatom=iatom+1 ica(i)=iatom - write (ipdb,10) iatom,restyp(iti),chainid(ichain),& + write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),& ires,(c(j,i),j=1,3) if (iti.ne.10) 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) 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) + if (itype(i,mnum).eq.ntyp1) cycle + if (itype(i,mnum).eq.10 .and. itype(i+1,mnum).ne.ntyp1_molec(mnum)) 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,mnum).ne.ntyp1_molec(mnum)) 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,mnum).eq.ntyp1_molec(mnum)) then write (ipdb,30) ica(i),ica(i)+1 endif enddo - if (itype(nct).ne.10) then + if (itype(nct,molnum(nct)).ne.10) then write (ipdb,30) ica(nct),ica(nct)+1 endif do i=1,nss diff --git a/source/wham/wham_calc.f90 b/source/wham/wham_calc.f90 index b2d16eb..05a07b4 100644 --- a/source/wham/wham_calc.f90 +++ b/source/wham/wham_calc.f90 @@ -342,6 +342,8 @@ 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 DEBUG write (iout,'(3i5,6f5.2,14f12.3)') i,ib,iparm,(ft(l),l=1,6),& evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,& @@ -379,7 +381,7 @@ +wturn3*eello_turn3 & +wturn6*eello_turn6+wel_loc*eel_loc & +edihcnstr+wtor_d*etors_d+wsccor*esccor & - +wbond*estr + +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation #else etot=wsc*evdw+wscp*evdw2 & +welec*(ees+evdw1) & @@ -389,7 +391,7 @@ +wturn3*eello_turn3 & +wturn6*eello_turn6+wel_loc*eel_loc+edihcnstr & +wtor_d*etors_d+wsccor*esccor & - +wbond*estr + +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation #endif #ifdef DEBUG @@ -752,6 +754,9 @@ ! edihcnstr=enetb(20,t,iparm) edihcnstr=enetb(19,t,iparm) edihcnstr=0.0d0 + ecationcation=enetb(42,t,iparm) + ecation_prot=enetb(41,t,iparm) + do k=0,nGridT betaT=startGridT+k*delta_T temper=betaT @@ -842,7 +847,7 @@ +ft(2)*wturn3*eello_turn3 & +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor & - +wbond*estr + +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees & +ftprim(1)*wtor*etors+ & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+ & @@ -865,7 +870,7 @@ +ft(2)*wturn3*eello_turn3 & +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor & - +wbond*estr + +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1) & +ftprim(1)*wtor*etors+ & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+ & diff --git a/source/wham/wham_data.f90 b/source/wham/wham_data.f90 index 1dfc3bc..f505680 100644 --- a/source/wham/wham_data.f90 +++ b/source/wham/wham_data.f90 @@ -1,7 +1,7 @@ module wham_data !--------------------------------------------------------------------------- !--------------------------------------------------------------------------- - integer,parameter :: max_eneW=21 + integer,parameter :: max_eneW=49 integer,parameter :: maxQ=1 integer,parameter :: maxQ1=MaxQ+2 integer,parameter :: max_parm=1 -- 1.7.9.5