X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fio_wham.F90;h=fee1e46c5912644b8b4a9966a5cf0c920860f4fc;hb=bc23440fbe68672d430f71f22f46b11265f003db;hp=97b74d122eba92ddcd66ce12386be27b8960a6ce;hpb=dacd7cb1d5f105094518f9d3b1428287fa9435f0;p=unres4.git diff --git a/source/wham/io_wham.F90 b/source/wham/io_wham.F90 index 97b74d1..fee1e46 100644 --- a/source/wham/io_wham.F90 +++ b/source/wham/io_wham.F90 @@ -98,6 +98,8 @@ call mygetenv('SCPPAR_NUCL',scpname_nucl) open (iscpp_nucl,file=scpname_nucl,status='old') + call mygetenv('IONPAR_NUCL',ionnuclname) + open (iionnucl,file=ionnuclname,status='old') #ifndef OLDSCP @@ -163,9 +165,9 @@ ! use energy_data use geometry_data, only:nres,deg2rad,c,dc,nres_molec - use control_data, only:iscode + use control_data, only:iscode,dyn_ss use io_base, only:rescode - use control, only:setup_var,init_int_table + use control, only:setup_var,init_int_table,hpb_partition use geometry, only:alloc_geo_arrays use energy, only:alloc_ener_arrays ! implicit real*8 (a-h,o-z) @@ -344,6 +346,20 @@ dhpb(i),ebr,forcon(i) enddo endif + if (ns.gt.0.and.dyn_ss) then + do i=nss+1,nhpb + ihpb(i-nss)=ihpb(i) + jhpb(i-nss)=jhpb(i) + forcon(i-nss)=forcon(i) + dhpb(i-nss)=dhpb(i) + enddo + nhpb=nhpb-nss + nss=0 + call hpb_partition + do i=1,ns + dyn_ss_mask(iss(i))=.true. + enddo + endif write (iout,'(a)') return end subroutine molread @@ -410,7 +426,7 @@ real(kind=8) :: dwa16,akl,si,rri,epsij,rrij,sigeps,sigt1sq,epsijlip,& sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm, & epspeptube,epssctube,sigmapeptube, & - sigmasctube + sigmasctube,ssscale real(kind=8) :: v0ij,v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,rjunk,& res1 @@ -434,6 +450,21 @@ call card_concat(controlcard,.true.) wname(4)="WCORRH" !el + call reada(controlcard,"D0CM",d0cm,3.78d0) + call reada(controlcard,"AKCM",akcm,15.1d0) + call reada(controlcard,"AKTH",akth,11.0d0) + call reada(controlcard,"AKCT",akct,12.0d0) + call reada(controlcard,"V1SS",v1ss,-1.08d0) + call reada(controlcard,"V2SS",v2ss,7.61d0) + call reada(controlcard,"V3SS",v3ss,13.7d0) + call reada(controlcard,"EBR",ebr,-5.50D0) + call reada(controlcard,"ATRISS",atriss,0.301D0) + call reada(controlcard,"BTRISS",btriss,0.021D0) + call reada(controlcard,"CTRISS",ctriss,1.001D0) + call reada(controlcard,"DTRISS",dtriss,1.001D0) + call reada(controlcard,"SSSCALE",ssscale,1.0D0) + dyn_ss=(index(controlcard,'DYN_SS').gt.0) + allocate(ww(max_eneW)) do i=1,n_eneW key = wname(i)(:ilen(wname(i))) @@ -485,6 +516,7 @@ allocate(ww(max_eneW)) wpepbase=ww(47) wscpho=ww(48) wpeppho=ww(49) + wcatnucl=ww(50) ! print *,"KURWA",ww(48) ! "WSCBASE ","WPEPBASE ","WSCPHO ","WPEPPHO " ! "WVDWPP ","WELPP ","WVDWPSB ","WELPSB ","WVDWSB ",& @@ -541,6 +573,8 @@ allocate(ww(max_eneW)) weights(47)= wpepbase weights(48) =wscpho weights(49) =wpeppho + weights(50) =wcatnucl + ! el-------- call card_concat(controlcard,.false.) @@ -584,7 +618,7 @@ allocate(ww(max_eneW)) write (iout,*) "Parameter set:",iparm write (iout,*) "Energy-term weights:" do i=1,n_eneW - write (iout,'(a16,f10.5)') wname(i),ww(i) + write (iout,'(i3,a16,f10.5)') i,wname(i),ww(i) enddo write (iout,*) "Sidechain potential file : ",& sidename_t(:ilen(sidename_t)) @@ -683,6 +717,23 @@ allocate(ww(max_eneW)) enddo print *, catprm endif + allocate(catnuclprm(14,ntyp_molec(2),ntyp_molec(5))) + do i=1,ntyp_molec(5) + do j=1,ntyp_molec(2) + write(iout,*) i,j + read(iionnucl,*) (catnuclprm(k,j,i),k=1,14) + enddo + enddo + write(*,'(3(5x,a6)11(7x,a6))') "w1 ","w2 ","epslj ","pis1 ", & + "sigma0","epsi0 ","chi1 ","chip1 ","sig ","b1 ","b2 ", & + "b3 ","b4 ","chis1 " + do i=1,ntyp_molec(5) + do j=1,ntyp_molec(2) + write(*,'(3(f10.3,x),11(f12.6,x),a3,2a)') (catnuclprm(k,j,i),k=1,14), & + restyp(i,5),"-",restyp(j,2) + enddo + enddo + read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2) do i=1,ntyp_molec(2) nbondterm_nucl(i)=1 @@ -2677,6 +2728,7 @@ allocate(ww(max_eneW)) ! Ions by Aga allocate(alphapolcat(ntyp,ntyp),epsheadcat(ntyp,ntyp),sig0headcat(ntyp,ntyp)) + allocate(alphapolcat2(ntyp,ntyp)) allocate(sigiso1cat(ntyp,ntyp),rborn1cat(ntyp,ntyp),rborn2cat(ntyp,ntyp),sigmap1cat(ntyp,ntyp)) allocate(sigmap2cat(ntyp,ntyp),sigiso2cat(ntyp,ntyp)) allocate(chis1cat(ntyp,ntyp),chis2cat(ntyp,ntyp),wquadcat(ntyp,ntyp),chipp1cat(ntyp,ntyp),chipp2cat(ntyp,ntyp)) @@ -2733,7 +2785,7 @@ allocate(ww(max_eneW)) ! rborncat(i,j),rborncat(j,i),& rborn1cat(i,j),rborn2cat(i,j),& (wqdipcat(k,i,j),k=1,2), & - alphapolcat(i,j),alphapolcat(j,i), & + alphapolcat(i,j),alphapolcat2(j,i), & (alphisocat(k,i,j),k=1,4),sigiso1cat(i,j),sigiso2cat(i,j),epsintabcat(i,j),debaykapcat(i,j) ! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i) END DO @@ -2913,7 +2965,7 @@ allocate(ww(max_eneW)) ! ! Define the constants of the disulfide bridge ! - ebr=-5.50D0 +! ebr=-5.50D0 ! ! Old arbitrary potential - commented out. ! @@ -2924,14 +2976,28 @@ allocate(ww(max_eneW)) ! energy surface of diethyl disulfide. ! A. Liwo and U. Kozlowska, 11/24/03 ! - D0CM = 3.78d0 - AKCM = 15.1d0 - AKTH = 11.0d0 - AKCT = 12.0d0 - V1SS =-1.08d0 - V2SS = 7.61d0 - V3SS = 13.7d0 +! D0CM = 3.78d0 +! AKCM = 15.1d0 +! AKTH = 11.0d0 +! AKCT = 12.0d0 +! V1SS =-1.08d0 +! V2SS = 7.61d0 +! V3SS = 13.7d0 +#ifndef CLUSTER + if (dyn_ss) then + ss_depth=(ebr/wsc-0.25*eps(1,1))*ssscale + Ht=(Ht/wsc-0.25*eps(1,1))*ssscale + akcm=akcm/wsc*ssscale + akth=akth/wsc*ssscale + akct=akct/wsc*ssscale + v1ss=v1ss/wsc*ssscale + v2ss=v2ss/wsc*ssscale + v3ss=v3ss/wsc*ssscale + else + ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain + endif +#endif if (lprint) then write (iout,'(/a)') "Disulfide bridge parameters:" write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr @@ -3013,7 +3079,7 @@ allocate(ww(max_eneW)) use control_data, only:indpdb,symetr,r_cut_ele,rlamb_ele,ions,& scelemode,TUBEmode,tor_mode,energy_dec - use energy_data, only:distchainmax,tubeR0,tubecenter + use energy_data, only:distchainmax,tubeR0,tubecenter,dyn_ss use geometry_data, only:boxxsize,boxysize,boxzsize,bordtubetop,& bordtubebot,tubebufthick,buftubebot,buftubetop ! implicit none @@ -3098,7 +3164,7 @@ allocate(ww(max_eneW)) energy_dec=(index(controlcard,'ENERGY_DEC').gt.0) call readi(controlcard,"SCELEMODE",scelemode,0) call readi(controlcard,"OLDION",oldion,0) - + dyn_ss=(index(controlcard,'DYN_SS').gt.0) print *,"SCELE",scelemode call readi(controlcard,'TORMODE',tor_mode,0) !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then @@ -3842,8 +3908,9 @@ allocate(ww(max_eneW)) !-------------------------------------------------------------------------------- subroutine pdboutW(ii,temp,efree,etot,entropy,rmsdev) - use geometry_data, only:nres,c + use geometry_data, only:nres,c,boxxsize,boxysize,boxzsize use energy_data, only:nss,nnt,nct,ihpb,jhpb,itype,molnum + use energy, only:boxshift ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'DIMENSIONS.ZSCOPT' @@ -3858,7 +3925,7 @@ allocate(ww(max_eneW)) 'D','E','F','G','H','I','J','K','L','M','N','O',& 'P','Q','R','S','V','W','X','Y','Z'/),shape(chainid)) integer,dimension(nres) :: ica !(maxres) - real(kind=8) :: temp,efree,etot,entropy,rmsdev + real(kind=8) :: temp,efree,etot,entropy,rmsdev,xj,yj,zj integer :: ii,i,j,iti,ires,iatom,ichain,mnum write(ipdb,'("REMARK CONF",i8," TEMPERATURE",f7.1," RMS",0pf7.2)')& ii,temp,rmsdev @@ -3882,9 +3949,17 @@ allocate(ww(max_eneW)) ires=ires+1 iatom=iatom+1 ica(i)=iatom + if (mnum.ne.5) then write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),& ires,(c(j,i),j=1,3) - if (iti.ne.10) then + else + xj=boxshift(c(1,i)-c(1,2),boxxsize) + yj=boxshift(c(2,i)-c(2,2),boxysize) + zj=boxshift(c(3,i)-c(3,2),boxzsize) + write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),& + ires,c(1,2)+xj,c(2,2)+yj,c(3,2)+zj + endif + if ((iti.ne.10).and.(mnum.ne.5)) then iatom=iatom+1 write (ipdb,20) iatom,restyp(iti,mnum),chainid(ichain),& ires,(c(j,nres+i),j=1,3)