X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fio_wham.F90;h=fee1e46c5912644b8b4a9966a5cf0c920860f4fc;hb=b2ef3fdfb9fd19710ea49377832d6c3791f04c11;hp=b6ee2006b8e1c2685e4ca654248200f7d9ab50f1;hpb=c27f7da3393f9e3a5b9d46fba0bbd520490a9fcc;p=unres4.git diff --git a/source/wham/io_wham.F90 b/source/wham/io_wham.F90 index b6ee200..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,14 +426,14 @@ 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 integer :: i,j,ichir1,ichir2,k,l,m,kk,ii,mm,junk,lll,ll,llll,n,jj integer :: nlobi,iblock,maxinter,iscprol,ncatprotparm character*3 string - + ! ! Body ! @@ -429,10 +445,26 @@ vbl=3.8D0 vblinv=1.0D0/vbl vblinv2=vblinv*vblinv + itime_mat=0 #ifndef CLUSTER 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))) @@ -461,6 +493,7 @@ allocate(ww(max_eneW)) wscloc=ww(12) wtor=ww(13) wtor_d=ww(14) + wstrain=ww(15) wvdwpp=ww(16) wbond=ww(18) wsccor=ww(19) @@ -483,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 ",& @@ -512,7 +546,7 @@ allocate(ww(max_eneW)) weights(12)=wscloc weights(13)=wtor weights(14)=wtor_d - weights(15)=0 !wstrain ! + weights(15)=wstrain !0 weights(16)=wvdwpp ! weights(17)=wbond weights(18)=0 !scal14 ! @@ -536,8 +570,11 @@ allocate(ww(max_eneW)) weights(41) =wcatcat weights(42) =wcatprot weights(46) =wscbase - weights(47) =wscpho - weights(48) =wpeppho + weights(47)= wpepbase + weights(48) =wscpho + weights(49) =wpeppho + weights(50) =wcatnucl + ! el-------- call card_concat(controlcard,.false.) @@ -581,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)) @@ -665,6 +702,7 @@ allocate(ww(max_eneW)) endif if (.not. allocated(msc)) allocate(msc(ntyp1,5)) if (.not. allocated(restok)) allocate(restok(ntyp1,5)) + if (oldion.eq.1) then do i=1,ntyp_molec(5) read(iion,*) msc(i,5),restok(i,5) @@ -678,6 +716,24 @@ allocate(ww(max_eneW)) read (iion,*) (catprm(i,k),i=1,ncatprotparm) 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 @@ -2669,6 +2725,126 @@ allocate(ww(max_eneW)) bad(i,j)=0.0D0 enddo enddo +! 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)) + allocate(epsintabcat(ntyp,ntyp)) + allocate(dtailcat(2,ntyp,ntyp)) + allocate(alphasurcat(4,ntyp,ntyp),alphisocat(4,ntyp,ntyp)) + allocate(wqdipcat(2,ntyp,ntyp)) + allocate(wstatecat(4,ntyp,ntyp)) + allocate(dheadcat(2,2,ntyp,ntyp)) + allocate(nstatecat(ntyp,ntyp)) + allocate(debaykapcat(ntyp,ntyp)) + + if (.not.allocated(epscat)) allocate (epscat(0:ntyp1,0:ntyp1)) + if (.not.allocated(sigmacat)) allocate(sigmacat(0:ntyp1,0:ntyp1)) +! if (.not.allocated(chicat)) allocate(chicat(ntyp1,ntyp1)) !(ntyp,ntyp) + if (.not.allocated(chi1cat)) allocate(chi1cat(ntyp1,ntyp1)) !(ntyp,ntyp) + if (.not.allocated(chi2cat)) allocate(chi2cat(ntyp1,ntyp1)) !(ntyp,ntyp) + + + allocate (ichargecat(ntyp_molec(5))) +! i to SC, j to jon, isideocat - nazwa pliku z ktorego czytam parametry + if (oldion.eq.0) then + if (.not.allocated(icharge)) then ! this mean you are oprating in old sc-sc mode + allocate(icharge(1:ntyp1)) + read(iion,*) (icharge(i),i=1,ntyp) + else + read(iion,*) ijunk + endif + + do i=1,ntyp_molec(5) + read(iion,*) msc(i,5),restok(i,5),ichargecat(i) + print *,msc(i,5),restok(i,5) + enddo + ip(5)=0.2 +!DIR$ NOUNROLL + do j=1,ntyp_molec(5) + do i=1,ntyp +! do j=1,ntyp_molec(5) +! write (*,*) "Im in ALAB", i, " ", j + read(iion,*) & + epscat(i,j),sigmacat(i,j), & +! chicat(i,j),chicat(j,i),chippcat(i,j),chippcat(j,i), & + chi1cat(i,j),chi2cat(i,j),chipp1cat(i,j),chipp2cat(i,j), & + + (alphasurcat(k,i,j),k=1,4),sigmap1cat(i,j),sigmap2cat(i,j),& +! chiscat(i,j),chiscat(j,i), & + chis1cat(i,j),chis2cat(i,j), & + + nstatecat(i,j),(wstatecat(k,i,j),k=1,4), & !5 w tej lini - 1 integer pierwszy + dheadcat(1,1,i,j),dheadcat(1,2,i,j),dheadcat(2,1,i,j),dheadcat(2,2,i,j),& + dtailcat(1,i,j),dtailcat(2,i,j), & + epsheadcat(i,j),sig0headcat(i,j), & +!wdipcat = w1 , w2 +! rborncat(i,j),rborncat(j,i),& + rborn1cat(i,j),rborn2cat(i,j),& + (wqdipcat(k,i,j),k=1,2), & + 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 + END DO + allocate(aa_aq_cat(-ntyp:ntyp,ntyp),bb_aq_cat(-ntyp:ntyp,ntyp)) + do i=1,ntyp + do j=1,ntyp_molec(5) + epsij=epscat(i,j) + rrij=sigmacat(i,j) + rrij=rrij**expon + sigeps=dsign(1.0D0,epsij) + epsij=dabs(epsij) + aa_aq_cat(i,j)=epsij*rrij*rrij + bb_aq_cat(i,j)=-sigeps*epsij*rrij + enddo + enddo + do i=1,ntyp + do j=1,ntyp_molec(5) + if (i.eq.10) then + write (iout,*) 'i= ', i, ' j= ', j + write (iout,*) 'epsi0= ', epscat(i,j) + write (iout,*) 'sigma0= ', sigmacat(i,j) + write (iout,*) 'chi1= ', chi1cat(i,j) + write (iout,*) 'chi1= ', chi2cat(i,j) + write (iout,*) 'chip1= ', chipp1cat(1,j) + write (iout,*) 'chip2= ', chipp2cat(1,j) + write (iout,*) 'alphasur1= ', alphasurcat(1,1,j) + write (iout,*) 'alphasur2= ', alphasurcat(2,1,j) + write (iout,*) 'alphasur3= ', alphasurcat(3,1,j) + write (iout,*) 'alphasur4= ', alphasurcat(4,1,j) + write (iout,*) 'sig1= ', sigmap1cat(1,j) + write (iout,*) 'sig2= ', sigmap2cat(1,j) + write (iout,*) 'chis1= ', chis1cat(1,j) + write (iout,*) 'chis1= ', chis2cat(1,j) + write (iout,*) 'nstatecat(i,j)= ', nstatecat(1,j) + write (iout,*) 'wstatecat(k,i,j),k=1= ',wstatecat(1,1,j) + write (iout,*) 'dhead= ', dheadcat(1,1,1,j) + write (iout,*) 'dhead2= ', dheadcat(1,2,1,j) + write (iout,*) 'a1= ', rborn1cat(i,j) + write (iout,*) 'a2= ', rborn2cat(i,j) + write (iout,*) 'epsin= ', epsintabcat(1,j), epsintabcat(j,1) + write (iout,*) 'alphapol1= ', alphapolcat(1,j) + write (iout,*) 'alphapol2= ', alphapolcat(j,1) + write (iout,*) 'w1= ', wqdipcat(1,i,j) + write (iout,*) 'w2= ', wqdipcat(2,i,j) + write (iout,*) 'debaykapcat(i,j)= ', debaykapcat(1,j) + endif + + If ((i.eq.1).and.(j.eq.27)) then + write (iout,*) 'SEP' + Write (iout,*) 'w1= ', wqdipcat(1,1,27) + Write (iout,*) 'w2= ', wqdipcat(2,1,27) + endif + + enddo + enddo + + endif + #ifdef CLUSTER ! ! Define the SC-p interaction constants @@ -2789,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. ! @@ -2800,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 @@ -2887,9 +3077,9 @@ allocate(ww(max_eneW)) subroutine read_general_data(*) use control_data, only:indpdb,symetr,r_cut_ele,rlamb_ele,ions,& - scelemode,TUBEmode,tor_mode + 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 @@ -2971,7 +3161,10 @@ allocate(ww(max_eneW)) call reada(controlcard,'BOXX',boxxsize,100.0d0) call reada(controlcard,'BOXY',boxysize,100.0d0) call reada(controlcard,'BOXZ',boxzsize,100.0d0) + 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 @@ -2991,7 +3184,7 @@ allocate(ww(max_eneW)) buftubetop=bordtubetop-tubebufthick endif ions=index(controlcard,"IONS").gt.0 - call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0) + call reada(controlcard,"R_CUT_ELE",r_cut_ele,25.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 @@ -3715,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' @@ -3731,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 @@ -3755,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)