X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fio_wham.F90;h=187502ee4a77211a6c83eecc930861497c007ddb;hb=5ec51dd877ddaf7847f1ce15c490656a063dfceb;hp=958e244e69d472adeb0b60d77b36997d80453129;hpb=db2040af09e73f3fc871158afe11ee242ec6cfb1;p=unres4.git diff --git a/source/wham/io_wham.F90 b/source/wham/io_wham.F90 index 958e244..187502e 100644 --- a/source/wham/io_wham.F90 +++ b/source/wham/io_wham.F90 @@ -51,6 +51,8 @@ ! Get parameter filenames and open the parameter files. call mygetenv('BONDPAR',bondname) open (ibond,file=bondname,status='old') + call mygetenv('BONDPAR_NUCL',bondname_nucl) + open (ibond_nucl,file=bondname_nucl,status='old') call mygetenv('THETPAR',thetname) open (ithep,file=thetname,status='old') call mygetenv('ROTPAR',rotname) @@ -69,6 +71,35 @@ open (isidep,file=sidename,status='old') call mygetenv('SIDEP',sidepname) open (isidep1,file=sidepname,status="old") + call mygetenv('THETPAR_NUCL',thetname_nucl) + open (ithep_nucl,file=thetname_nucl,status='old') + call mygetenv('ROTPAR_NUCL',rotname_nucl) + open (irotam_nucl,file=rotname_nucl,status='old') + call mygetenv('TORPAR_NUCL',torname_nucl) + open (itorp_nucl,file=torname_nucl,status='old') + call mygetenv('TORDPAR_NUCL',tordname_nucl) + open (itordp_nucl,file=tordname_nucl,status='old') + call mygetenv('SIDEPAR_NUCL',sidename_nucl) + open (isidep_nucl,file=sidename_nucl,status='old') + call mygetenv('SIDEPAR_SCBASE',sidename_scbase) + open (isidep_scbase,file=sidename_scbase,status='old') + call mygetenv('PEPPAR_PEPBASE',pepname_pepbase) + open (isidep_pepbase,file=pepname_pepbase,status='old') + call mygetenv('SCPAR_PHOSPH',pepname_scpho) + open (isidep_scpho,file=pepname_scpho,status='old') + call mygetenv('PEPPAR_PHOSPH',pepname_peppho) + open (isidep_peppho,file=pepname_peppho,status='old') + call mygetenv('LIPTRANPAR',liptranname) + open (iliptranpar,file=liptranname,status='old',action='read') + call mygetenv('TUBEPAR',tubename) + open (itube,file=tubename,status='old',action='read') + call mygetenv('IONPAR',ionname) + open (iion,file=ionname,status='old',action='read') + + call mygetenv('SCPPAR_NUCL',scpname_nucl) + open (iscpp_nucl,file=scpname_nucl,status='old') + + #ifndef OLDSCP ! ! 8/9/01 In the newest version SCp interaction constants are read from a file @@ -377,7 +408,10 @@ integer :: iparm,nkcctyp !el real(kind=8) :: ip,mp real(kind=8) :: dwa16,akl,si,rri,epsij,rrij,sigeps,sigt1sq,epsijlip,& - sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm + sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm, & + epspeptube,epssctube,sigmapeptube, & + sigmasctube + 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 @@ -427,11 +461,39 @@ 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) wcatcat=ww(42) wcatprot=ww(41) + wcorr3_nucl=ww(38) + wcorr_nucl=ww(37) + wtor_d_nucl=ww(36) + wtor_nucl=ww(35) + wsbloc=ww(34) + wang_nucl=ww(33) + wbond_nucl=ww(32) + welsb=ww(31) + wvdwsb=ww(30) + welpsb=ww(29) + wvdwpsb=ww(28) + welpp=ww(27) + wvdwpp_nucl=ww(26) + wscbase=ww(46) + wpepbase=ww(47) + wscpho=ww(48) + wpeppho=ww(49) +! print *,"KURWA",ww(48) +! "WSCBASE ","WPEPBASE ","WSCPHO ","WPEPPHO " +! "WVDWPP ","WELPP ","WVDWPSB ","WELPSB ","WVDWSB ",& +! "WELSB ","WBOND ","WANG ","WSBLOC ","WTOR ",& +! "WTORD ","WCORR ","WCORR3 ","WNULL ","WNULL ",& +! "WCATPROT ","WCATCAT +! +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl& +! +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb& +! +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl& +! +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl& endif ! @@ -451,14 +513,32 @@ allocate(ww(max_eneW)) weights(12)=wscloc weights(13)=wtor weights(14)=wtor_d - weights(15)=0 !wstrain ! - weights(16)=0 !wvdwpp ! + weights(15)=wstrain !0 + weights(16)=wvdwpp ! weights(17)=wbond weights(18)=0 !scal14 ! weights(21)=wsccor weights(42)=wcatprot weights(41)=wcatcat - + weights(26)= wvdwpp_nucl + + weights(27) =welpp + weights(28) =wvdwpsb + weights(29) =welpsb + weights(30) =wvdwsb + weights(31) =welsb + weights(32) =wbond_nucl + weights(33) =wang_nucl + weights(34) =wsbloc + weights(35) =wtor_nucl + weights(36) =wtor_d_nucl + weights(37) =wcorr_nucl + weights(38) =wcorr3_nucl + weights(41) =wcatcat + weights(42) =wcatprot + weights(46) =wscbase + weights(47) =wscpho + weights(48) =wpeppho ! el-------- call card_concat(controlcard,.false.) @@ -536,6 +616,10 @@ allocate(ww(max_eneW)) allocate(nbondterm(ntyp)) !(ntyp) allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp) allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp) + allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp) + allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp) + allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp) + !el allocate(msc(ntyp+1)) !(ntyp+1) !el allocate(isc(ntyp+1)) !(ntyp+1) !el allocate(restok(ntyp+1)) !(ntyp+1) @@ -595,6 +679,17 @@ allocate(ww(max_eneW)) read (iion,*) (catprm(i,k),i=1,ncatprotparm) enddo print *, catprm + read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2) + do i=1,ntyp_molec(2) + nbondterm_nucl(i)=1 + read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),rjunk,restok(i,2) +! dsc(i) = vbldsc0_nucl(1,i) +! if (i.eq.10) then +! dsc_inv(i)=0.0D0 +! else +! dsc_inv(i)=1.0D0/dsc(i) +! endif + enddo !---------------------------------------------------- @@ -933,6 +1028,67 @@ allocate(ww(max_eneW)) ENDIF ! TOR_MODE #endif +!--------------- Reading theta parameters for nucleic acid------- + read (ithep_nucl,*) nthetyp_nucl,ntheterm_nucl,& + ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl + nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl) + allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1) + allocate(aa0thet_nucl(nthetyp_nucl+1,& + nthetyp_nucl+1,nthetyp_nucl+1)) +!(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2) + allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,& + nthetyp_nucl+1,nthetyp_nucl+1)) +!(maxtheterm,-maxthetyp1:maxthetyp1,& +! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2) + allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,& + nthetyp_nucl+1,nthetyp_nucl+1)) + allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,& + nthetyp_nucl+1,nthetyp_nucl+1)) + allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,& + nthetyp_nucl+1,nthetyp_nucl+1)) + allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,& + nthetyp_nucl+1,nthetyp_nucl+1)) +!(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,& +! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2) + allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,& + nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1)) + allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,& + nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1)) + +!(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,& +! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)) + + read (ithep_nucl,*) (ithetyp_nucl(i),i=1,ntyp1_molec(2)) + + aa0thet_nucl(:,:,:)=0.0d0 + aathet_nucl(:,:,:,:)=0.0d0 + bbthet_nucl(:,:,:,:,:)=0.0d0 + ccthet_nucl(:,:,:,:,:)=0.0d0 + ddthet_nucl(:,:,:,:,:)=0.0d0 + eethet_nucl(:,:,:,:,:)=0.0d0 + ffthet_nucl(:,:,:,:,:,:)=0.0d0 + ggthet_nucl(:,:,:,:,:,:)=0.0d0 + + do i=1,nthetyp_nucl + do j=1,nthetyp_nucl + do k=1,nthetyp_nucl + read (ithep_nucl,'(3a)') t1,t2,t3 + read (ithep_nucl,*) aa0thet_nucl(i,j,k) + read (ithep_nucl,*)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl) + read (ithep_nucl,*) & + (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), & + (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), & + (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), & + (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl) + read (ithep_nucl,*) & + (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), & + ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), & + llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl) + enddo + enddo + enddo + + !------------------------------------------- allocate(nlob(ntyp1)) !(ntyp1) allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp) @@ -1070,6 +1226,24 @@ allocate(ww(max_eneW)) enddo #endif close(irotam) +!---------reading nucleic acid parameters for rotamers------------------- + allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp) + do i=1,ntyp_molec(2) + read (irotam_nucl,*) + do j=1,9 + read(irotam_nucl,*) sc_parmin_nucl(j,i) + enddo + enddo + close(irotam_nucl) + if (lprint) then + write (iout,*) + write (iout,*) "Base rotamer parameters" + do i=1,ntyp_molec(2) + write (iout,'(a)') restyp(i,2) + write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9) + enddo + endif + read (ifourier,*) nloctyp !el write(iout,*)"nloctyp",nloctyp @@ -1651,6 +1825,11 @@ allocate(ww(max_eneW)) enddo enddo endif +#ifndef NEWCORR + do i=1,ntyp1 + itype2loc(i)=itortyp(i) + enddo +#endif ELSE IF (TOR_MODE.eq.1) THEN !C read valence-torsional parameters @@ -2117,6 +2296,10 @@ allocate(ww(max_eneW)) allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp) if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1)) allocate(r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1) + allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),& + dcavtub(ntyp1)) + allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),& + tubetranene(ntyp1)) do i=1,ntyp1 do j=1,ntyp1 aa_aq(i,j)=0.0D0 @@ -2215,6 +2398,271 @@ allocate(ww(max_eneW)) endif enddo enddo + allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2))) + allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2)) + allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2)) + allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2) + allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2)) + allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2) + allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2)) + allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2) + allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2) + allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2) + allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2) + allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2) + allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2) + allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2) + allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2)) + allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2) + +! augm(:,:)=0.0D0 +! chip(:)=0.0D0 +! alp(:)=0.0D0 +! sigma0(:)=0.0D0 +! sigii(:)=0.0D0 +! rr0(:)=0.0D0 + + read (isidep_nucl,*) ipot_nucl +! print *,"TU?!",ipot_nucl + if (ipot_nucl.eq.1) then + do i=1,ntyp_molec(2) + do j=i,ntyp_molec(2) + read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),& + elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j) + enddo + enddo + else + do i=1,ntyp_molec(2) + do j=i,ntyp_molec(2) + read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),& + chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),& + elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j) + enddo + enddo + endif +! rpp(1,1)=2**(1.0/6.0)*5.16158 + do i=1,ntyp_molec(2) + do j=i,ntyp_molec(2) + rrij=sigma_nucl(i,j) + r0_nucl(i,j)=rrij + r0_nucl(j,i)=rrij + rrij=rrij**expon + epsij=4*eps_nucl(i,j) + sigeps=dsign(1.0D0,epsij) + epsij=dabs(epsij) + aa_nucl(i,j)=epsij*rrij*rrij + bb_nucl(i,j)=-sigeps*epsij*rrij + ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij) + ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij + ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij + ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij + sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- & + 2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i))) + enddo + do j=1,i-1 + aa_nucl(i,j)=aa_nucl(j,i) + bb_nucl(i,j)=bb_nucl(j,i) + ael3_nucl(i,j)=ael3_nucl(j,i) + ael6_nucl(i,j)=ael6_nucl(j,i) + ael63_nucl(i,j)=ael63_nucl(j,i) + ael32_nucl(i,j)=ael32_nucl(j,i) + elpp3_nucl(i,j)=elpp3_nucl(j,i) + elpp6_nucl(i,j)=elpp6_nucl(j,i) + elpp63_nucl(i,j)=elpp63_nucl(j,i) + elpp32_nucl(i,j)=elpp32_nucl(j,i) + eps_nucl(i,j)=eps_nucl(j,i) + sigma_nucl(i,j)=sigma_nucl(j,i) + sigmaii_nucl(i,j)=sigmaii_nucl(j,i) + enddo + enddo + + write(iout,*) "tube param" + read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, & + ccavtubpep,dcavtubpep,tubetranenepep + sigmapeptube=sigmapeptube**6 + sigeps=dsign(1.0D0,epspeptube) + epspeptube=dabs(epspeptube) + pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2 + pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube + write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep + do i=1,ntyp + read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), & + ccavtub(i),dcavtub(i),tubetranene(i) + sigmasctube=sigmasctube**6 + sigeps=dsign(1.0D0,epssctube) + epssctube=dabs(epssctube) + sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2 + sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube + write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i) + enddo +!-----------------READING SC BASE POTENTIALS----------------------------- + allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2))) + allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2))) + allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2)) + allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2)) + allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2))) + allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2))) + allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2))) + allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2)) + allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2))) + allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2))) + allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2))) + allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2))) + allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2))) + allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2))) + allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2))) + allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2))) + + + do i=1,ntyp_molec(1) + do j=1,ntyp_molec(2)-1 ! without U then we will take T for U + write (*,*) "Im in ", i, " ", j + read(isidep_scbase,*) & + eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),& + chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2) + write(*,*) "eps",eps_scbase(i,j) + read(isidep_scbase,*) & + (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), & + chis_scbase(i,j,1),chis_scbase(i,j,2) + read(isidep_scbase,*) & + dhead_scbasei(i,j), & + dhead_scbasej(i,j), & + rborn_scbasei(i,j),rborn_scbasej(i,j) + read(isidep_scbase,*) & + (wdipdip_scbase(k,i,j),k=1,3), & + (wqdip_scbase(k,i,j),k=1,2) + read(isidep_scbase,*) & + alphapol_scbase(i,j), & + epsintab_scbase(i,j) + END DO + END DO + allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2))) + allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2))) + + do i=1,ntyp_molec(1) + do j=1,ntyp_molec(2)-1 + epsij=eps_scbase(i,j) + rrij=sigma_scbase(i,j) +! r0(i,j)=rrij +! r0(j,i)=rrij + rrij=rrij**expon +! epsij=eps(i,j) + sigeps=dsign(1.0D0,epsij) + epsij=dabs(epsij) + aa_scbase(i,j)=epsij*rrij*rrij + bb_scbase(i,j)=-sigeps*epsij*rrij + enddo + enddo +!-----------------READING PEP BASE POTENTIALS------------------- + allocate(eps_pepbase(ntyp_molec(2))) + allocate(sigma_pepbase(ntyp_molec(2))) + allocate(chi_pepbase(ntyp_molec(2),2)) + allocate(chipp_pepbase(ntyp_molec(2),2)) + allocate(alphasur_pepbase(4,ntyp_molec(2))) + allocate(sigmap1_pepbase(ntyp_molec(2))) + allocate(sigmap2_pepbase(ntyp_molec(2))) + allocate(chis_pepbase(ntyp_molec(2),2)) + allocate(wdipdip_pepbase(3,ntyp_molec(2))) + + + do j=1,ntyp_molec(2)-1 ! without U then we will take T for U + write (*,*) "Im in ", i, " ", j + read(isidep_pepbase,*) & + eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),& + chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2) + write(*,*) "eps",eps_pepbase(j) + read(isidep_pepbase,*) & + (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), & + chis_pepbase(j,1),chis_pepbase(j,2) + read(isidep_pepbase,*) & + (wdipdip_pepbase(k,j),k=1,3) + END DO + allocate(aa_pepbase(ntyp_molec(2))) + allocate(bb_pepbase(ntyp_molec(2))) + + do j=1,ntyp_molec(2)-1 + epsij=eps_pepbase(j) + rrij=sigma_pepbase(j) +! r0(i,j)=rrij +! r0(j,i)=rrij + rrij=rrij**expon +! epsij=eps(i,j) + sigeps=dsign(1.0D0,epsij) + epsij=dabs(epsij) + aa_pepbase(j)=epsij*rrij*rrij + bb_pepbase(j)=-sigeps*epsij*rrij + enddo +!--------------READING SC PHOSPHATE------------------------------------- +!--------------READING SC PHOSPHATE------------------------------------- + allocate(eps_scpho(ntyp_molec(1))) + allocate(sigma_scpho(ntyp_molec(1))) + allocate(chi_scpho(ntyp_molec(1),2)) + allocate(chipp_scpho(ntyp_molec(1),2)) + allocate(alphasur_scpho(4,ntyp_molec(1))) + allocate(sigmap1_scpho(ntyp_molec(1))) + allocate(sigmap2_scpho(ntyp_molec(1))) + allocate(chis_scpho(ntyp_molec(1),2)) + allocate(wqq_scpho(ntyp_molec(1))) + allocate(wqdip_scpho(2,ntyp_molec(1))) + allocate(alphapol_scpho(ntyp_molec(1))) + allocate(epsintab_scpho(ntyp_molec(1))) + allocate(dhead_scphoi(ntyp_molec(1))) + allocate(rborn_scphoi(ntyp_molec(1))) + allocate(rborn_scphoj(ntyp_molec(1))) + allocate(alphi_scpho(ntyp_molec(1))) + + +! j=1 + do j=1,ntyp_molec(1) ! without U then we will take T for U + write (*,*) "Im in scpho ", i, " ", j + read(isidep_scpho,*) & + eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),& + chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2) + write(*,*) "eps",eps_scpho(j) + read(isidep_scpho,*) & + (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), & + chis_scpho(j,1),chis_scpho(j,2) + read(isidep_scpho,*) & + (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j) + read(isidep_scpho,*) & + epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), & + alphi_scpho(j) + + END DO + allocate(aa_scpho(ntyp_molec(1))) + allocate(bb_scpho(ntyp_molec(1))) + do j=1,ntyp_molec(1) + epsij=eps_scpho(j) + rrij=sigma_scpho(j) +! r0(i,j)=rrij +! r0(j,i)=rrij + rrij=rrij**expon +! epsij=eps(i,j) + sigeps=dsign(1.0D0,epsij) + epsij=dabs(epsij) + aa_scpho(j)=epsij*rrij*rrij + bb_scpho(j)=-sigeps*epsij*rrij + enddo + + + read(isidep_peppho,*) & + eps_peppho,sigma_peppho + read(isidep_peppho,*) & + (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho + read(isidep_peppho,*) & + (wqdip_peppho(k),k=1,2) + + epsij=eps_peppho + rrij=sigma_peppho +! r0(i,j)=rrij +! r0(j,i)=rrij + rrij=rrij**expon +! epsij=eps(i,j) + sigeps=dsign(1.0D0,epsij) + epsij=dabs(epsij) + aa_peppho=epsij*rrij*rrij + bb_peppho=-sigeps*epsij*rrij + allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2) do i=1,ntyp @@ -2233,6 +2681,48 @@ allocate(ww(max_eneW)) enddo enddo #endif + allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1) + read (itorp_nucl,*) ntortyp_nucl +! print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2) +!el from energy module--------- + allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2) + allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2) + + allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor) + allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) + allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor) + allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2) + + allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) + allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) +!el--------------------------- + nterm_nucl(:,:)=0 + nlor_nucl(:,:)=0 +!el-------------------- + read (itorp_nucl,*) & + (itortyp_nucl(i),i=1,ntyp_molec(2)) +! print *,itortyp_nucl(:) +!c write (iout,*) 'ntortyp',ntortyp + do i=1,ntortyp_nucl + do j=1,ntortyp_nucl + read (itorp_nucl,*) nterm_nucl(i,j),nlor_nucl(i,j) +! print *,nterm_nucl(i,j),nlor_nucl(i,j) + v0ij=0.0d0 + si=-1.0d0 + do k=1,nterm_nucl(i,j) + read (itorp_nucl,*) kk,v1_nucl(k,i,j),v2_nucl(k,i,j) + v0ij=v0ij+si*v1_nucl(k,i,j) + si=-si + enddo + do k=1,nlor_nucl(i,j) + read (itorp_nucl,*) kk,vlor1_nucl(k,i,j),& + vlor2_nucl(k,i,j),vlor3_nucl(k,i,j) + v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2) + enddo + v0_nucl(i,j)=v0ij + enddo + enddo + !elwrite(iout,*) "parmread kontrol before oldscp" ! @@ -2283,6 +2773,20 @@ allocate(ww(max_eneW)) enddo endif #endif + allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2) + + do i=1,ntyp_molec(2) + read (iscpp_nucl,*) eps_scp_nucl(i),rscp_nucl(i) + enddo + do i=1,ntyp_molec(2) + aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12 + bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6 + enddo + r0pp=1.12246204830937298142*5.16158 + epspp=4.95713/4 + AEES=108.661 + BEES=0.433246 + ! ! Define the constants of the disulfide bridge ! @@ -3109,6 +3613,7 @@ allocate(ww(max_eneW)) do k=1,3 cref(k,j+i,kkk)=cref(k,j,kkk) enddo + write(iout,*) "J",j,"J+I",j+i phi_ref(j+i)=phi_ref(j) theta_ref(j+i)=theta_ref(j) alph_ref(j+i)=alph_ref(j) @@ -3223,8 +3728,9 @@ allocate(ww(max_eneW)) ! include 'COMMON.HEADER' ! include 'COMMON.SBRIDGE' character(len=50) :: tytul - character(len=1),dimension(10) :: chainid=reshape((/'A','B','C',& - 'D','E','F','G','H','I','J'/),shape(chainid)) + character(len=1),dimension(24) :: chainid=reshape((/'A','B','C',& + '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 integer :: ii,i,j,iti,ires,iatom,ichain,mnum @@ -3241,9 +3747,11 @@ allocate(ww(max_eneW)) mnum=molnum(i) iti=itype(i,mnum) if (iti.eq.ntyp1) then + if (itype(i-1,molnum(i-1)).eq.ntyp1) then ichain=ichain+1 ires=0 write (ipdb,'(a)') 'TER' + endif else ires=ires+1 iatom=iatom+1