X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;ds=sidebyside;f=source%2Fwham%2Fio_wham.F90;h=d8cb743260e13037c45ce9283b1d1ad709d79742;hb=ae2f9918e1e5f88a78cc32b1d3a7ff90dede1207;hp=4c8a65db35f0ca711f4c1862b52cb2e089c8dda7;hpb=b74333aa5a167dd84f48d8a9a4c7ac62058e9132;p=unres4.git diff --git a/source/wham/io_wham.F90 b/source/wham/io_wham.F90 index 4c8a65d..d8cb743 100644 --- a/source/wham/io_wham.F90 +++ b/source/wham/io_wham.F90 @@ -100,6 +100,8 @@ open (iscpp_nucl,file=scpname_nucl,status='old') call mygetenv('IONPAR_NUCL',ionnuclname) open (iionnucl,file=ionnuclname,status='old') + call mygetenv('IONPAR_TRAN',iontranname) + open (iiontran,file=iontranname,status='old',action='read') #ifndef OLDSCP @@ -490,7 +492,7 @@ call reada(controlcard,"DTRISS",dtriss,1.001D0) call reada(controlcard,"SSSCALE",ssscale,1.0D0) dyn_ss=(index(controlcard,'DYN_SS').gt.0) - + call reada(controlcard,"LIPSCALE",lipscale,1.0D0) allocate(ww(max_eneW)) do i=1,n_eneW key = wname(i)(:ilen(wname(i))) @@ -543,6 +545,7 @@ allocate(ww(max_eneW)) wscpho=ww(48) wpeppho=ww(49) wcatnucl=ww(50) + wcat_tran=ww(56) ! print *,"KURWA",ww(48) ! "WSCBASE ","WPEPBASE ","WSCPHO ","WPEPPHO " ! "WVDWPP ","WELPP ","WVDWPSB ","WELPSB ","WVDWSB ",& @@ -600,6 +603,7 @@ allocate(ww(max_eneW)) weights(48) =wscpho weights(49) =wpeppho weights(50) =wcatnucl + weights(56)=wcat_tran ! el-------- call card_concat(controlcard,.false.) @@ -726,8 +730,8 @@ allocate(ww(max_eneW)) enddo enddo endif - if (.not. allocated(msc)) allocate(msc(ntyp1,5)) - if (.not. allocated(restok)) allocate(restok(ntyp1,5)) + if (.not. allocated(msc)) allocate(msc(-ntyp1:ntyp1,5)) + if (.not. allocated(restok)) allocate(restok(-ntyp1:ntyp1,5)) if (oldion.eq.1) then do i=1,ntyp_molec(5) @@ -2774,7 +2778,7 @@ allocate(ww(max_eneW)) if (.not.allocated(chi2cat)) allocate(chi2cat(ntyp1,ntyp1)) !(ntyp,ntyp) - allocate (ichargecat(ntyp_molec(5))) + allocate (ichargecat(-ntyp_molec(5):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 @@ -2784,13 +2788,13 @@ allocate(ww(max_eneW)) read(iion,*) ijunk endif - do i=1,ntyp_molec(5) + do i=-ntyp_molec(5),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 j=1,ntyp_molec(5)-1 do i=1,ntyp ! do j=1,ntyp_molec(5) ! write (*,*) "Im in ALAB", i, " ", j @@ -2818,7 +2822,7 @@ allocate(ww(max_eneW)) 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) + do j=1,ntyp_molec(5)-1 !without zinc epsij=epscat(i,j) rrij=sigmacat(i,j) rrij=rrij**expon @@ -2870,7 +2874,36 @@ allocate(ww(max_eneW)) enddo endif - +! here two denotes the Zn2+ and Cu2+ + write(iout,*) "before TRANPARM" + allocate(aomicattr(0:3,2)) + allocate(athetacattran(0:6,5,2)) + allocate(agamacattran(3,5,2)) + allocate(acatshiftdsc(5,2)) + allocate(bcatshiftdsc(5,2)) + allocate(demorsecat(5,2)) + allocate(alphamorsecat(5,2)) + allocate(x0catleft(5,2)) + allocate(x0catright(5,2)) + allocate(x0cattrans(5,2)) + allocate(ntrantyp(2)) + do i=1,1 ! currently only Zn2+ + + read(iiontran,*) ntrantyp(i) +!ntrantyp=4 +!| ao0 ao1 ao2 ao3 +!ASP| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -.5 +!CYS| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0left x0right x0transi +!GLU| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -0.5 +!HIS| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -.5 + read(iiontran,*) (aomicattr(j,i),j=0,3) + do j=1,ntrantyp(i) + read (iiontran,*) (agamacattran(k,j,i),k=1,3),& + (athetacattran(k,j,i),k=0,6),acatshiftdsc(j,i),bcatshiftdsc(j,i),& + demorsecat(j,i),alphamorsecat(j,i),x0catleft(j,i),x0catright(j,i),& + x0cattrans(j,i) + enddo + enddo #ifdef CLUSTER ! ! Define the SC-p interaction constants @@ -3103,11 +3136,13 @@ allocate(ww(max_eneW)) subroutine read_general_data(*) use control_data, only:indpdb,symetr,r_cut_ele,rlamb_ele,ions,& - scelemode,TUBEmode,tor_mode,energy_dec + scelemode,TUBEmode,tor_mode,energy_dec,r_cut_ang,r_cut_mart,& + rlamb_mart use energy_data, only:distchainmax,tubeR0,tubecenter,dyn_ss,constr_homology use geometry_data, only:boxxsize,boxysize,boxzsize,bordtubetop,& - bordtubebot,tubebufthick,buftubebot,buftubetop + bordtubebot,tubebufthick,buftubebot,buftubetop,buflipbot, bufliptop,bordlipbot,bordliptop, & + lipbufthick,lipthick ! implicit none ! include "DIMENSIONS" ! include "DIMENSIONS.ZSCOPT" @@ -3126,6 +3161,7 @@ allocate(ww(max_eneW)) ! include "COMMON.FREE" ! include "COMMON.CONTROL" ! include "COMMON.ENERGIES" + character(len=800) :: controlcard integer :: i,j,k,ii,n_ene_found integer :: ind,itype1,itype2,itypf,itypsc,itypp @@ -3187,6 +3223,23 @@ 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) + call reada(controlcard,"LIPTHICK",lipthick,0.0d0) + call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0) + if (lipthick.gt.0.0d0) then + bordliptop=(boxzsize+lipthick)/2.0 + bordlipbot=bordliptop-lipthick + if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) & + write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE" + buflipbot=bordlipbot+lipbufthick + bufliptop=bordliptop-lipbufthick + if ((lipbufthick*2.0d0).gt.lipthick) & + write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF" + endif !lipthick.gt.0 + write(iout,*) "bordliptop=",bordliptop + write(iout,*) "bordlipbot=",bordlipbot + write(iout,*) "bufliptop=",bufliptop + write(iout,*) "buflipbot=",buflipbot + energy_dec=(index(controlcard,'ENERGY_DEC').gt.0) call readi(controlcard,"SCELEMODE",scelemode,0) call readi(controlcard,"OLDION",oldion,0) @@ -3220,6 +3273,9 @@ allocate(ww(max_eneW)) 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 + call reada(controlcard,"R_CUT_MART",r_cut_mart,15.0d0) + call reada(controlcard,"LAMBDA_MART",rlamb_mart,0.3d0) + call reada(controlcard,"R_CUT_ANG",r_cut_ang,4.2d0) check_conf=index(controlcard,"NO_CHECK_CONF").eq.0 call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0) call readi(controlcard,'SYM',symetr,1) @@ -3972,8 +4028,8 @@ allocate(ww(max_eneW)) do i=nnt,nct mnum=molnum(i) iti=itype(i,mnum) - if (iti.eq.ntyp1) then - if (itype(i-1,molnum(i-1)).eq.ntyp1) then + if (iti.eq.ntyp1_molec(mnum)) then + if (itype(i-1,molnum(i-1)).eq.ntyp1_molec(mnum)) then ichain=ichain+1 ires=0 write (ipdb,'(a)') 'TER' @@ -4002,7 +4058,7 @@ allocate(ww(max_eneW)) write (ipdb,'(a)') 'TER' do i=nnt,nct-1 mnum=molnum(i) - if (itype(i,mnum).eq.ntyp1) cycle + if (itype(i,mnum).eq.ntyp1_molec(mnum)) 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,mnum).ne.10 .and. itype(i+1,mnum).ne.ntyp1_molec(mnum)) then