X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fio_wham.f90;h=eab5a49c8560a7d9ee6273f6350edad1e6e0686e;hb=705644e0cbb7678faefd6fe1bc436159d38ad85d;hp=399dc1a2ac278c24d9d607547c9f995950242384;hpb=3d82fd4865ef986a636dcb144357cc095a20ffbc;p=unres4.git 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