X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fio_config.f90;h=87ddf94beea67828febd369586e764f506f5e616;hb=10689ab7d813dfbdbb0c6e631d90234d78ea306a;hp=0dce2e11063e1dd91b449b27860c8455841cf329;hpb=5d299c1a16ab51f8206b8ee3b17c7bcabe9321b7;p=unres4.git diff --git a/source/unres/io_config.f90 b/source/unres/io_config.f90 index 0dce2e1..87ddf94 100644 --- a/source/unres/io_config.f90 +++ b/source/unres/io_config.f90 @@ -478,7 +478,7 @@ ! write (iout,100) ! do i=1,nres -! write (iout,110) restyp(itype(i)),i,c(1,i),c(2,i),& +! write (iout,110) restyp(itype(i,1)),i,c(1,i),c(2,i),& ! c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i) ! enddo ! 100 format (//' alpha-carbon coordinates ',& @@ -529,9 +529,9 @@ ! read secondary structure prediction from JPRED here! ! read(isecpred,'(A80)',err=100,end=100) line ! read(line,'(f10.3)',err=110) ftors - read(isecpred,'(f10.3)',err=110) ftors + read(isecpred,'(f10.3)',err=110) ftors(1) - write (iout,*) 'FTORS factor =',ftors + write (iout,*) 'FTORS factor =',ftors(1) ! initialize secstruc to any do i=1,nres secstruc(i) ='-' @@ -553,6 +553,7 @@ ii=0 do i=1,nres + ftors(i)=ftors(1) if ( secstruc(i) .eq. 'H') then ! Helix restraints for this residue ii=ii+1 @@ -735,15 +736,16 @@ character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/) logical :: lprint,LaTeX real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob) - real(kind=8),dimension(13) :: b + real(kind=8),dimension(13) :: buse character(len=3) :: lancuch !,ucase !el local variables integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm - integer :: maxinter,junk,kk,ii + integer :: maxinter,junk,kk,ii,ncatprotparm real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,& dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,& sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,& - res1,epsijlip + res1,epsijlip,epspeptube,epssctube,sigmapeptube, & + sigmasctube integer :: ichir1,ichir2 ! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) !el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)) @@ -775,12 +777,12 @@ ! allocate(dsc(ntyp1)) !(ntyp1) allocate(dsc_inv(ntyp1)) !(ntyp1) + 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) allocate(nbondterm(ntyp)) !(ntyp) allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp) allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp) - allocate(msc(ntyp+1)) !(ntyp+1) - allocate(isc(ntyp+1)) !(ntyp+1) - allocate(restok(ntyp+1)) !(ntyp+1) allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp) allocate(long_r_sidechain(ntyp)) allocate(short_r_sidechain(ntyp)) @@ -788,6 +790,10 @@ dsc_inv(:)=0.0d0 #ifdef CRYST_BOND + allocate(msc(ntyp+1)) !(ntyp+1) + allocate(isc(ntyp+1)) !(ntyp+1) + allocate(restok(ntyp+1)) !(ntyp+1) + read (ibond,*) vbldp0,akp,mp,ip,pstok do i=1,ntyp nbondterm(i)=1 @@ -800,10 +806,19 @@ endif enddo #else - read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp,ip,pstok - do i=1,ntyp + mp(:)=0.0d0 + ip(:)=0.0d0 + msc(:,:)=0.0d0 + isc(:,:)=0.0d0 + + allocate(msc(ntyp+1,5)) !(ntyp+1) + allocate(isc(ntyp+1,5)) !(ntyp+1) + allocate(restok(ntyp+1,5)) !(ntyp+1) + + read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp(1),ip(1),pstok(1) + do i=1,ntyp_molec(1) read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),& - j=1,nbondterm(i)),msc(i),isc(i),restok(i) + j=1,nbondterm(i)),msc(i,1),isc(i,1),restok(i,1) dsc(i) = vbldsc0(1,i) if (i.eq.10) then dsc_inv(i)=0.0D0 @@ -812,20 +827,57 @@ endif enddo #endif + 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),isc(i,2),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 +! read (ibond_nucl,*) junk,vbldp0_nucl,akp_nucl,rjunk,mp(2),ip(2),pstok(2) +! do i=1,ntyp_molec(2) +! read (ibond_nucl,*) nbondterm_nucl(i),(vbldsc0_nucl(j,i),& +! aksc_nucl(j,i),abond0_nucl(j,i),& +! j=1,nbondterm_nucl(i)),msc(i,2),isc(i,2),restok(i,2) +! dsc(i) = vbldsc0(1,i) +! if (i.eq.10) then +! dsc_inv(i)=0.0D0 +! else +! dsc_inv(i)=1.0D0/dsc(i) +! endif +! enddo + if (lprint) then write(iout,'(/a/)')"Dynamic constants of the interaction sites:" write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',& 'inertia','Pstok' - write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok + write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp(1),ip(1),pstok(1) do i=1,ntyp - write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),& - vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i) + write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),& + vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i,1),isc(i,1),restok(i,1) do j=2,nbondterm(i) write (iout,'(13x,3f10.5)') & vbldsc0(j,i),aksc(j,i),abond0(j,i) enddo enddo endif + 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 +! read (iion,*) (vcatprm(k),k=1,ncatprotpram) !---------------------------------------------------- allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp)) allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp) @@ -911,7 +963,7 @@ ' ATHETA0 ',' A1 ',' A2 ',& ' B1 ',' B2 ' do i=1,ntyp - write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,& + write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,& a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2) enddo write (iout,'(/a/9x,5a/79(1h-))') & @@ -919,7 +971,7 @@ ' ALPH0 ',' ALPH1 ',' ALPH2 ',& ' ALPH3 ',' SIGMA0C ' do i=1,ntyp - write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,& + write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,& (polthet(j,i),j=0,3),sigc0(i) enddo write (iout,'(/a/9x,5a/79(1h-))') & @@ -927,7 +979,7 @@ ' THETA0 ',' SIGMA0 ',' G1 ',& ' G2 ',' G3 ' do i=1,ntyp - write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),& + write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,theta0(i),& sig0(i),(gthet(j,i),j=1,3) enddo else @@ -938,7 +990,7 @@ ' theta0 ',' a1*10^2 ',' a2*10^2 ',& ' b1*10^1 ',' b2*10^1 ' do i=1,ntyp - write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),& + write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),& a0thet(i),(100*athet(j,i,1,1),j=1,2),& (10*bthet(j,i,1,1),j=1,2) enddo @@ -947,7 +999,7 @@ ' alpha0 ',' alph1 ',' alph2 ',& ' alhp3 ',' sigma0c ' do i=1,ntyp - write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),& + write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i,1),& (polthet(j,i),j=0,3),sigc0(i) enddo write (iout,'(/a/9x,5a/79(1h-))') & @@ -955,7 +1007,7 @@ ' theta0 ',' sigma0*10^2 ',' G1*10^-1',& ' G2 ',' G3*10^1 ' do i=1,ntyp - write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),& + write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),theta0(i),& 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0 enddo endif @@ -1206,6 +1258,65 @@ close (ithep_pdb) #endif close(ithep) +!--------------- Reading theta parameters for nucleic acid------- + read (ithep_nucl,*,err=111,end=111) 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,*,err=111,end=111) (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)',end=111,err=111) t1,t2,t3 + read (ithep_nucl,*,end=111,err=111) aa0thet_nucl(i,j,k) + read (ithep_nucl,*,end=111,err=111)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl) + read (ithep_nucl,*,end=111,err=111) & + (((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,*,end=111,err=111) & + (((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) @@ -1284,7 +1395,7 @@ nlobi=nlob(i) if (nlobi.gt.0) then if (LaTeX) then - write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),& + write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i,1),& ' # of gaussian lobes:',nlobi,' dsc:',dsc(i) write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') & 'log h',(bsc(j,i),j=1,nlobi) @@ -1325,6 +1436,24 @@ enddo endif enddo +!---------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,*,end=112,err=112) + do j=1,9 + read(irotam_nucl,*,end=112,err=112) 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 the parameters of the probability distribution/energy expression ! of the side chains. @@ -1594,6 +1723,48 @@ enddo endif #endif + allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1) + read (itorp_nucl,*,end=113,err=113) 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,*,end=113,err=113) & + (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,*,end=113,err=113) 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,*,end=113,err=113) 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,*,end=113,err=113) 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 + ! Read of Side-chain backbone correlation parameters ! Modified 11 May 2012 by Adasko !CC @@ -1737,6 +1908,7 @@ #endif if (lprint) then + l=3 write (iout,'(/a/)') 'Torsional constants:' do i=1,nsccortyp do j=1,nsccortyp @@ -1781,88 +1953,88 @@ do i=0,nloctyp-1 read (ifourier,*,end=115,err=115) - read (ifourier,*,end=115,err=115) (b(ii),ii=1,13) + read (ifourier,*,end=115,err=115) (buse(ii),ii=1,13) if (lprint) then write (iout,*) 'Type',i - write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13) + write (iout,'(a,i2,a,f10.5)') ('buse(',ii,')=',buse(ii),ii=1,13) endif - B1(1,i) = b(3) - B1(2,i) = b(5) - B1(1,-i) = b(3) - B1(2,-i) = -b(5) -! b1(1,i)=0.0d0 -! b1(2,i)=0.0d0 - B1tilde(1,i) = b(3) - B1tilde(2,i) =-b(5) - B1tilde(1,-i) =-b(3) - B1tilde(2,-i) =b(5) -! b1tilde(1,i)=0.0d0 -! b1tilde(2,i)=0.0d0 - B2(1,i) = b(2) - B2(2,i) = b(4) - B2(1,-i) =b(2) - B2(2,-i) =-b(4) - -! b2(1,i)=0.0d0 -! b2(2,i)=0.0d0 - CC(1,1,i)= b(7) - CC(2,2,i)=-b(7) - CC(2,1,i)= b(9) - CC(1,2,i)= b(9) - CC(1,1,-i)= b(7) - CC(2,2,-i)=-b(7) - CC(2,1,-i)=-b(9) - CC(1,2,-i)=-b(9) + B1(1,i) = buse(3) + B1(2,i) = buse(5) + B1(1,-i) = buse(3) + B1(2,-i) = -buse(5) +! buse1(1,i)=0.0d0 +! buse1(2,i)=0.0d0 + B1tilde(1,i) = buse(3) + B1tilde(2,i) =-buse(5) + B1tilde(1,-i) =-buse(3) + B1tilde(2,-i) =buse(5) +! buse1tilde(1,i)=0.0d0 +! buse1tilde(2,i)=0.0d0 + B2(1,i) = buse(2) + B2(2,i) = buse(4) + B2(1,-i) =buse(2) + B2(2,-i) =-buse(4) + +! buse2(1,i)=0.0d0 +! buse2(2,i)=0.0d0 + CC(1,1,i)= buse(7) + CC(2,2,i)=-buse(7) + CC(2,1,i)= buse(9) + CC(1,2,i)= buse(9) + CC(1,1,-i)= buse(7) + CC(2,2,-i)=-buse(7) + CC(2,1,-i)=-buse(9) + CC(1,2,-i)=-buse(9) ! CC(1,1,i)=0.0d0 ! CC(2,2,i)=0.0d0 ! CC(2,1,i)=0.0d0 ! CC(1,2,i)=0.0d0 - Ctilde(1,1,i)=b(7) - Ctilde(1,2,i)=b(9) - Ctilde(2,1,i)=-b(9) - Ctilde(2,2,i)=b(7) - Ctilde(1,1,-i)=b(7) - Ctilde(1,2,-i)=-b(9) - Ctilde(2,1,-i)=b(9) - Ctilde(2,2,-i)=b(7) + Ctilde(1,1,i)=buse(7) + Ctilde(1,2,i)=buse(9) + Ctilde(2,1,i)=-buse(9) + Ctilde(2,2,i)=buse(7) + Ctilde(1,1,-i)=buse(7) + Ctilde(1,2,-i)=-buse(9) + Ctilde(2,1,-i)=buse(9) + Ctilde(2,2,-i)=buse(7) ! Ctilde(1,1,i)=0.0d0 ! Ctilde(1,2,i)=0.0d0 ! Ctilde(2,1,i)=0.0d0 ! Ctilde(2,2,i)=0.0d0 - DD(1,1,i)= b(6) - DD(2,2,i)=-b(6) - DD(2,1,i)= b(8) - DD(1,2,i)= b(8) - DD(1,1,-i)= b(6) - DD(2,2,-i)=-b(6) - DD(2,1,-i)=-b(8) - DD(1,2,-i)=-b(8) + DD(1,1,i)= buse(6) + DD(2,2,i)=-buse(6) + DD(2,1,i)= buse(8) + DD(1,2,i)= buse(8) + DD(1,1,-i)= buse(6) + DD(2,2,-i)=-buse(6) + DD(2,1,-i)=-buse(8) + DD(1,2,-i)=-buse(8) ! DD(1,1,i)=0.0d0 ! DD(2,2,i)=0.0d0 ! DD(2,1,i)=0.0d0 ! DD(1,2,i)=0.0d0 - Dtilde(1,1,i)=b(6) - Dtilde(1,2,i)=b(8) - Dtilde(2,1,i)=-b(8) - Dtilde(2,2,i)=b(6) - Dtilde(1,1,-i)=b(6) - Dtilde(1,2,-i)=-b(8) - Dtilde(2,1,-i)=b(8) - Dtilde(2,2,-i)=b(6) + Dtilde(1,1,i)=buse(6) + Dtilde(1,2,i)=buse(8) + Dtilde(2,1,i)=-buse(8) + Dtilde(2,2,i)=buse(6) + Dtilde(1,1,-i)=buse(6) + Dtilde(1,2,-i)=-buse(8) + Dtilde(2,1,-i)=buse(8) + Dtilde(2,2,-i)=buse(6) ! Dtilde(1,1,i)=0.0d0 ! Dtilde(1,2,i)=0.0d0 ! Dtilde(2,1,i)=0.0d0 ! Dtilde(2,2,i)=0.0d0 - EE(1,1,i)= b(10)+b(11) - EE(2,2,i)=-b(10)+b(11) - EE(2,1,i)= b(12)-b(13) - EE(1,2,i)= b(12)+b(13) - EE(1,1,-i)= b(10)+b(11) - EE(2,2,-i)=-b(10)+b(11) - EE(2,1,-i)=-b(12)+b(13) - EE(1,2,-i)=-b(12)-b(13) + EE(1,1,i)= buse(10)+buse(11) + EE(2,2,i)=-buse(10)+buse(11) + EE(2,1,i)= buse(12)-buse(13) + EE(1,2,i)= buse(12)+buse(13) + EE(1,1,-i)= buse(10)+buse(11) + EE(2,2,-i)=-buse(10)+buse(11) + EE(2,1,-i)=-buse(12)+buse(13) + EE(1,2,-i)=-buse(12)-buse(13) ! ee(1,1,i)=1.0d0 ! ee(2,2,i)=1.0d0 @@ -1926,6 +2098,7 @@ allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp) allocate(augm(ntyp,ntyp)) !(ntyp,ntyp) 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)) @@ -1964,7 +2137,7 @@ 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,1),sigma0(i),i=1,ntyp) endif ! goto 50 !----------------------- LJK potential -------------------------------- @@ -1978,13 +2151,15 @@ 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,1),sigma0(i),rr0(i),& i=1,ntyp) endif ! goto 50 !---------------------- GB or BP potential ----------------------------- case(3:4) ! 30 do i=1,ntyp +! print *,"I AM in SCELE",scelemode + if (scelemode.eq.0) then do i=1,ntyp read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp) enddo @@ -2009,9 +2184,84 @@ 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,1),sigma0(i),sigii(i),& chip(i),alp(i),i=1,ntyp) endif + else +! print *,ntyp,"NTYP" + allocate(icharge(ntyp1)) +! print *,ntyp,icharge(i) + icharge(:)=0 + read (isidep,*) (icharge(i),i=1,ntyp) + print *,ntyp,icharge(i) +! if(.not.allocated(eps)) allocate(eps(-ntyp +!c write (2,*) "icharge",(icharge(i),i=1,ntyp) + allocate(alphapol(ntyp,ntyp),epshead(ntyp,ntyp),sig0head(ntyp,ntyp)) + allocate(sigiso1(ntyp,ntyp),rborn(ntyp,ntyp),sigmap1(ntyp,ntyp)) + allocate(sigmap2(ntyp,ntyp),sigiso2(ntyp,ntyp)) + allocate(chis(ntyp,ntyp),wquad(ntyp,ntyp),chipp(ntyp,ntyp)) + allocate(epsintab(ntyp,ntyp)) + allocate(dtail(2,ntyp,ntyp)) + allocate(alphasur(4,ntyp,ntyp),alphiso(4,ntyp,ntyp)) + allocate(wqdip(2,ntyp,ntyp)) + allocate(wstate(4,ntyp,ntyp)) + allocate(dhead(2,2,ntyp,ntyp)) + allocate(nstate(ntyp,ntyp)) + if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1)) + if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp) + do i=1,ntyp + do j=1,i +! write (*,*) "Im in ALAB", i, " ", j + read(isidep,*) & + eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & + (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), & + chis(i,j),chis(j,i), & + nstate(i,j),(wstate(k,i,j),k=1,4), & + dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),& + dtail(1,i,j),dtail(2,i,j), & + epshead(i,j),sig0head(i,j), & + rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), & + alphapol(i,j),alphapol(j,i), & + (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j) +! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i) + END DO + END DO + DO i = 1, ntyp + DO j = i+1, ntyp + eps(i,j) = eps(j,i) + sigma(i,j) = sigma(j,i) + sigmap1(i,j)=sigmap1(j,i) + sigmap2(i,j)=sigmap2(j,i) + sigiso1(i,j)=sigiso1(j,i) + sigiso2(i,j)=sigiso2(j,i) +! print *,"ATU",sigma(j,i),sigma(i,j),i,j + nstate(i,j) = nstate(j,i) + dtail(1,i,j) = dtail(1,j,i) + dtail(2,i,j) = dtail(2,j,i) + DO k = 1, 4 + alphasur(k,i,j) = alphasur(k,j,i) + wstate(k,i,j) = wstate(k,j,i) + alphiso(k,i,j) = alphiso(k,j,i) + END DO + + dhead(2,1,i,j) = dhead(1,1,j,i) + dhead(2,2,i,j) = dhead(1,2,j,i) + dhead(1,1,i,j) = dhead(2,1,j,i) + dhead(1,2,i,j) = dhead(2,2,j,i) + + epshead(i,j) = epshead(j,i) + sig0head(i,j) = sig0head(j,i) + + DO k = 1, 2 + wqdip(k,i,j) = wqdip(k,j,i) + END DO + + wquad(i,j) = wquad(j,i) + epsintab(i,j) = epsintab(j,i) +! if (epsintab(i,j).ne.1.0) print *,"WHAT?",i,j,epsintab(i,j) + END DO + END DO + endif ! goto 50 !--------------------- GBV potential ----------------------------------- case(5) @@ -2026,7 +2276,7 @@ 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,1),sigma0(i),rr0(i),& sigii(i),chip(i),alp(i),i=1,ntyp) endif case default @@ -2035,21 +2285,36 @@ end select continue close (isidep) + !----------------------------------------------------------------------- ! Calculate the "working" parameters of SC interactions. !el from module energy - COMMON.INTERACT------- - allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp) + allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1)) + if (.not.allocated(chi)) allocate(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) + 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)) aa_aq(:,:)=0.0D0 bb_aq(:,:)=0.0D0 aa_lip(:,:)=0.0D0 bb_lip(:,:)=0.0D0 + if (scelemode.eq.0) then chi(:,:)=0.0D0 sigma(:,:)=0.0D0 r0(:,:)=0.0D0 - + endif + acavtub(:)=0.0d0 + bcavtub(:)=0.0d0 + ccavtub(:)=0.0d0 + dcavtub(:)=0.0d0 + sc_aa_tube_par(:)=0.0d0 + sc_bb_tube_par(:)=0.0d0 + !-------------------------------- do i=2,ntyp @@ -2058,6 +2323,7 @@ epslip(i,j)=epslip(j,i) enddo enddo + if (scelemode.eq.0) then do i=1,ntyp do j=i,ntyp sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2) @@ -2066,6 +2332,7 @@ rs0(j,i)=rs0(i,j) enddo enddo + endif if (lprint) write (iout,'(/a/10x,7a/72(1h-))') & 'Working parameters of the SC interactions:',& ' a ',' b ',' augm ',' sigma ',' r0 ',& @@ -2075,6 +2342,7 @@ epsij=eps(i,j) if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then rrij=sigma(i,j) +! print *,"SIGMA ZLA?",sigma(i,j) else rrij=rr0(i)+rr0(j) endif @@ -2085,6 +2353,7 @@ sigeps=dsign(1.0D0,epsij) epsij=dabs(epsij) aa_aq(i,j)=epsij*rrij*rrij + print *,"ADASKO",epsij,rrij,expon bb_aq(i,j)=-sigeps*epsij*rrij aa_aq(j,i)=aa_aq(i,j) bb_aq(j,i)=bb_aq(i,j) @@ -2096,7 +2365,7 @@ aa_lip(j,i)=aa_lip(i,j) bb_lip(j,i)=bb_lip(i,j) !C write(iout,*) aa_lip - if (ipot.gt.2) then + if ((ipot.gt.2).and. (scelemode.eq.0)) then sigt1sq=sigma0(i)**2 sigt2sq=sigma0(j)**2 sigii1=sigii(i) @@ -2128,12 +2397,278 @@ endif if (lprint) then write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') & - restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),& + restyp(i,1),restyp(j,1),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 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------------------------------------- + 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) bad(:,:)=0.0D0 @@ -2186,6 +2721,20 @@ endif ! lprint=.false. #endif + allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2) + + do i=1,ntyp_molec(2) + read (iscpp_nucl,*,end=118,err=118) 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 ! @@ -2222,6 +2771,31 @@ write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,& ' v3ss:',v3ss endif + if (shield_mode.gt.0) then + pi=4.0D0*datan(1.0D0) +!C VSolvSphere the volume of solving sphere + print *,pi,"pi" +!C rpp(1,1) is the energy r0 for peptide group contact and will be used for it +!C there will be no distinction between proline peptide group and normal peptide +!C group in case of shielding parameters + VSolvSphere=4.0/3.0*pi*(4.50d0)**3 + VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3 + write (iout,*) VSolvSphere,VSolvSphere_div +!C long axis of side chain + do i=1,ntyp + long_r_sidechain(i)=vbldsc0(1,i) +! if (scelemode.eq.0) then + short_r_sidechain(i)=sigma(i,i)/sqrt(2.0) + if (short_r_sidechain(i).eq.0.0) short_r_sidechain(i)=0.2 +! else +! short_r_sidechain(i)=sigma(i,i) +! endif + write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),& + sigma0(i) + enddo + buff_shield=1.0d0 + endif + return 111 write (iout,*) "Error reading bending energy parameters." goto 999 @@ -2280,11 +2854,11 @@ ! Read the PDB file and convert the peptide geometry into virtual-chain ! geometry. use geometry_data - use energy_data, only: itype + use energy_data, only: itype,istype use control_data use compare_data use MPI_data - use control, only: rescode +! use control, only: rescode,sugarcode ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.LOCAL' @@ -2297,7 +2871,7 @@ ! include 'COMMON.CONTROL' ! include 'COMMON.DISTFIT' ! include 'COMMON.SETUP' - integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift!,ity!,& + integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,& ! ishift_pdb logical :: lprn=.true.,fail real(kind=8),dimension(3) :: e1,e2,e3 @@ -2306,32 +2880,51 @@ character(len=5) :: atom character(len=80) :: card real(kind=8),dimension(3,20) :: sccor - integer :: kkk,lll,icha,kupa !rescode, + integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode, + integer :: isugar,molecprev,firstion + character*1 :: sugar real(kind=8) :: cou + real(kind=8),dimension(3) :: ccc !el local varables integer,dimension(2,maxres/3) :: hfrag_alloc integer,dimension(4,maxres/3) :: bfrag_alloc real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm) - + real(kind=8),dimension(:,:), allocatable :: c_temporary + integer,dimension(:,:) , allocatable :: itype_temporary + integer,dimension(:),allocatable :: istype_temp efree_temp=0.0d0 ibeg=1 ishift1=0 ishift=0 + molecule=1 + counter=0 ! write (2,*) "UNRES_PDB",unres_pdb ires=0 ires_old=0 +#ifdef WHAM_RUN + do i=1,nres + do j=1,5 + itype(i,j)=0 + enddo + enddo +#endif nres=0 iii=0 lsecondary=.false. nhfrag=0 nbfrag=0 + do j=1,5 + nres_molec(j)=0 + enddo + + !----------------------------- allocate(hfrag(2,maxres/3)) !(2,maxres/3) allocate(bfrag(4,maxres/3)) !(4,maxres/3) - + if(.not. allocated(istype)) allocate(istype(maxres)) do i=1,100000 read (ipdbin,'(a80)',end=10) card -! write (iout,'(a)') card + write (iout,'(a)') card if (card(:5).eq.'HELIX') then nhfrag=nhfrag+1 lsecondary=.true. @@ -2354,9 +2947,11 @@ else if (card(:3).eq.'TER') then ! End current chain ires_old=ires+2 + ishift=ishift+1 ishift1=ishift1+1 - itype(ires_old)=ntyp1 - itype(ires_old-1)=ntyp1 + itype(ires_old,molecule)=ntyp1_molec(molecule) + itype(ires_old-1,molecule)=ntyp1_molec(molecule) + nres_molec(molecule)=nres_molec(molecule)+2 ibeg=2 ! write (iout,*) "Chain ended",ires,ishift,ires_old if (unres_pdb) then @@ -2372,6 +2967,7 @@ ! Read free energy if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp ! Fish out the ATOM cards. +! write(iout,*) 'card',card(1:20) if (index(card(1:4),'ATOM').gt.0) then read (card(12:16),*) atom ! write (iout,*) "! ",atom," !",ires @@ -2405,7 +3001,8 @@ ishift=ires-1 if (res.ne.'GLY' .and. res.ne. 'ACE') then ishift=ishift-1 - itype(1)=ntyp1 + itype(1,1)=ntyp1 + nres_molec(molecule)=nres_molec(molecule)+1 endif ires=ires-ishift+ishift1 ires_old=ires @@ -2418,17 +3015,37 @@ ishift1=ishift1-1 !!!!! ! write (iout,*) "New chain started",ires,ishift,ishift1,"!" ires=ires-ishift+ishift1 +! print *,ires,ishift,ishift1 ires_old=ires ibeg=0 else ishift=ishift-(ires-ishift+ishift1-ires_old-1) ires=ires-ishift+ishift1 ires_old=ires - endif + endif +! print *,'atom',ires,atom if (res.eq.'ACE' .or. res.eq.'NHE') then - itype(ires)=10 + itype(ires,1)=10 + else + if (atom.eq.'CA '.or.atom.eq.'N ') then + molecule=1 + itype(ires,molecule)=rescode(ires,res,0,molecule) + firstion=0 +! nres_molec(molecule)=nres_molec(molecule)+1 else - itype(ires)=rescode(ires,res,0) + molecule=2 + itype(ires,molecule)=rescode(ires,res(2:3),0,molecule) +! nres_molec(molecule)=nres_molec(molecule)+1 + read (card(19:19),'(a1)') sugar + isugar=sugarcode(sugar,ires) +! if (ibeg.eq.1) then +! istype(1)=isugar +! else + istype(ires)=isugar +! print *,"ires=",ires,istype(ires) +! endif + + endif endif else ires=ires-ishift+ishift1 @@ -2441,91 +3058,197 @@ if (atom.eq.'CA' .or. atom.eq.'CH3' .or. & res.eq.'NHE'.and.atom(:2).eq.'HN') then read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3) +! print *,ires,ishift,ishift1 ! write (iout,*) "backbone ",atom #ifdef DEBUG write (iout,'(2i3,2x,a,3f8.3)') & - ires,itype(ires),res,(c(j,ires),j=1,3) + ires,itype(ires,1),res,(c(j,ires),j=1,3) #endif iii=iii+1 + nres_molec(molecule)=nres_molec(molecule)+1 do j=1,3 sccor(j,iii)=c(j,ires) enddo -! write (*,*) card(23:27),ires,itype(ires) + else if (.not.unres_pdb .and. (atom.eq."C1'" .or. & + atom.eq."C2'" .or. atom.eq."C3'" & + .or. atom.eq."C4'" .or. atom.eq."O4'")) then + read(card(31:54),'(3f8.3)') (ccc(j),j=1,3) +!c write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3) +! print *,ires,ishift,ishift1 + counter=counter+1 +! iii=iii+1 +! do j=1,3 +! sccor(j,iii)=c(j,ires) +! enddo + do j=1,3 + c(j,ires)=c(j,ires)+ccc(j)/5.0 + enddo + print *,counter,molecule + if (counter.eq.5) then +! iii=iii+1 + nres_molec(molecule)=nres_molec(molecule)+1 + firstion=0 +! do j=1,3 +! sccor(j,iii)=c(j,ires) +! enddo + counter=0 + endif +! print *, "ATOM",atom(1:3) + else if (atom.eq."C5'") then + read (card(19:19),'(a1)') sugar + isugar=sugarcode(sugar,ires) + if (ibeg.eq.1) then + istype(1)=isugar + else + istype(ires)=isugar +! print *,ires,istype(ires) + endif + if (unres_pdb) then + molecule=2 +! print *,"nres_molec(molecule)",nres_molec(molecule),ires + read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3) + nres_molec(molecule)=nres_molec(molecule)+1 + print *,"nres_molec(molecule)",nres_molec(molecule),ires + + else + iii=iii+1 + read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3) + endif + else if ((atom.eq."C1'").and.unres_pdb) then + iii=iii+1 + read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3) +! write (*,*) card(23:27),ires,itype(ires,1) else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. & atom.ne.'N' .and. atom.ne.'C' .and. & atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. & - atom.ne.'OXT' .and. atom(:2).ne.'3H') then + atom.ne.'OXT' .and. atom(:2).ne.'3H' & + .and. atom.ne.'P '.and. & + atom(1:1).ne.'H' .and. & + atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'& + ) then ! write (iout,*) "sidechain ",atom +! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3) + if ((molecule.ne.2).or.(atom(3:3).ne."'")) then +! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3) iii=iii+1 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3) + endif endif - endif + else if ((ions).and.(card(1:6).eq.'HETATM')) then + if (firstion.eq.0) then + firstion=1 + if (unres_pdb) then + do j=1,3 + dc(j,ires)=sccor(j,iii) + enddo + else + call sccenter(ires,iii,sccor) + endif + endif + read (card(12:16),*) atom + print *,"HETATOM", atom + read (card(18:20),'(a3)') res + if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.& + (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') & + .or.(atom(1:2).eq.'K ')) & + then + ires=ires+1 + if (molecule.ne.5) molecprev=molecule + molecule=5 + nres_molec(molecule)=nres_molec(molecule)+1 + print *,"HERE",nres_molec(molecule) + res=res(2:3)//' ' + itype(ires,molecule)=rescode(ires,res,0,molecule) + read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3) + endif + endif !atom enddo 10 write (iout,'(a,i5)') ' Number of residues found: ',ires if (ires.eq.0) return ! Calculate dummy residue coordinates inside the "chain" of a multichain ! system nres=ires + if (((ires_old.ne.ires).and.(molecule.ne.5)) & + ) & + nres_molec(molecule)=nres_molec(molecule)-2 + print *,'I have',nres, nres_molec(:) + + do k=1,4 ! ions are without dummy + if (nres_molec(k).eq.0) cycle do i=2,nres-1 -! write (iout,*) i,itype(i) -! if (itype(i).eq.ntyp1) then -! write (iout,*) "dummy",i,itype(i) +! write (iout,*) i,itype(i,1) +! if (itype(i,1).eq.ntyp1) then +! write (iout,*) "dummy",i,itype(i,1) ! do j=1,3 ! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2 ! dc(j,i)=c(j,i) ! enddo ! endif - if (itype(i).eq.ntyp1) then - if (itype(i+1).eq.ntyp1) then + if (itype(i,k).eq.ntyp1_molec(k)) then + if (itype(i+1,k).eq.ntyp1_molec(k)) then + if (itype(i+2,k).eq.0) then +! print *,"masz sieczke" + do j=1,5 + if (itype(i+2,j).ne.0) then + itype(i+1,k)=0 + itype(i+1,j)=ntyp1_molec(j) + nres_molec(k)=nres_molec(k)-1 + nres_molec(j)=nres_molec(j)+1 + go to 3331 + endif + enddo + 3331 continue + endif ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain -! first is connected prevous chain (itype(i+1).eq.ntyp1)=true -! second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false - if (unres_pdb) then +! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true +! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false +! if (unres_pdb) then ! 2/15/2013 by Adam: corrected insertion of the last dummy residue ! print *,i,'tu dochodze' - call refsys(i-3,i-2,i-1,e1,e2,e3,fail) - if (fail) then - e2(1)=0.0d0 - e2(2)=1.0d0 - e2(3)=0.0d0 - endif !fail - print *,i,'a tu?' - do j=1,3 - c(j,i)=c(j,i-1)-1.9d0*e2(j) - enddo - else !unres_pdb +! call refsys(i-3,i-2,i-1,e1,e2,e3,fail) +! if (fail) then +! e2(1)=0.0d0 +! e2(2)=1.0d0 +! e2(3)=0.0d0 +! endif !fail +! print *,i,'a tu?' +! do j=1,3 +! c(j,i)=c(j,i-1)-1.9d0*e2(j) +! enddo +! else !unres_pdb do j=1,3 dcj=(c(j,i-2)-c(j,i-3))/2.0 if (dcj.eq.0) dcj=1.23591524223 c(j,i)=c(j,i-1)+dcj c(j,nres+i)=c(j,i) enddo - endif !unres_pdb - else !itype(i+1).eq.ntyp1 - if (unres_pdb) then +! endif !unres_pdb + else !itype(i+1,1).eq.ntyp1 +! if (unres_pdb) then ! 2/15/2013 by Adam: corrected insertion of the first dummy residue - call refsys(i+1,i+2,i+3,e1,e2,e3,fail) - if (fail) then - e2(1)=0.0d0 - e2(2)=1.0d0 - e2(3)=0.0d0 - endif +! call refsys(i+1,i+2,i+3,e1,e2,e3,fail) +! if (fail) then +! e2(1)=0.0d0 +! e2(2)=1.0d0 +! e2(3)=0.0d0 +! endif do j=1,3 c(j,i)=c(j,i+1)-1.9d0*e2(j) enddo - else !unres_pdb +! else !unres_pdb do j=1,3 dcj=(c(j,i+3)-c(j,i+2))/2.0 if (dcj.eq.0) dcj=1.23591524223 c(j,i)=c(j,i+1)-dcj c(j,nres+i)=c(j,i) enddo - endif !unres_pdb - endif !itype(i+1).eq.ntyp1 +! endif !unres_pdb + endif !itype(i+1,1).eq.ntyp1 endif !itype.eq.ntyp1 enddo + enddo ! Calculate the CM of the last side chain. if (iii.gt.0) then if (unres_pdb) then @@ -2539,28 +3262,33 @@ ! nres=ires nsup=nres nstart_sup=1 - if (itype(nres).ne.10) then +! print *,"molecule",molecule + if ((itype(nres,1).ne.10)) then nres=nres+1 - itype(nres)=ntyp1 - if (unres_pdb) then + if (molecule.eq.5) molecule=molecprev + itype(nres,molecule)=ntyp1_molec(molecule) + nres_molec(molecule)=nres_molec(molecule)+1 +! if (unres_pdb) then ! 2/15/2013 by Adam: corrected insertion of the last dummy residue - call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail) - if (fail) then - e2(1)=0.0d0 - e2(2)=1.0d0 - e2(3)=0.0d0 - endif - do j=1,3 - c(j,nres)=c(j,nres-1)-1.9d0*e2(j) - enddo - else +! call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail) +! if (fail) then +! e2(1)=0.0d0 +! e2(2)=1.0d0 +! e2(3)=0.0d0 +! endif +! do j=1,3 +! c(j,nres)=c(j,nres-1)-1.9d0*e2(j) +! enddo +! else do j=1,3 dcj=(c(j,nres-2)-c(j,nres-3))/2.0 c(j,nres)=c(j,nres-1)+dcj c(j,2*nres)=c(j,nres) enddo - endif +! endif endif +! print *,'I have',nres, nres_molec(:) + !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb #ifdef WHAM_RUN if (nres.ne.nres0) then @@ -2608,7 +3336,8 @@ c(j,nres+1)=c(j,1) c(j,2*nres)=c(j,nres) enddo - if (itype(1).eq.ntyp1) then + + if (itype(1,1).eq.ntyp1) then nsup=nsup-1 nstart_sup=2 if (unres_pdb) then @@ -2630,6 +3359,24 @@ enddo endif endif +! First lets assign correct dummy to correct type of chain +! 1) First residue + if (itype(1,1).eq.ntyp1) then + if (itype(2,1).eq.0) then + do j=2,5 + if (itype(2,j).ne.0) then + itype(1,1)=0 + itype(1,j)=ntyp1_molec(j) + nres_molec(1)=nres_molec(1)-1 + nres_molec(j)=nres_molec(j)+1 + go to 3231 + endif + enddo +3231 continue + endif + endif + print *,'I have',nres, nres_molec(:) + ! Copy the coordinates to reference coordinates ! do i=1,2*nres ! do j=1,3 @@ -2640,11 +3387,11 @@ if (lprn) then write (iout,'(/a)') & "Cartesian coordinates of the reference structure" - write (iout,'(a,3(3x,a5),5x,3(3x,a5))') & + write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)" do ires=1,nres - write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') & - restyp(itype(ires)),ires,(c(j,ires),j=1,3),& + write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') & + (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),& (c(j,ires+nres),j=1,3) enddo endif @@ -2654,12 +3401,83 @@ write (iout,'(a)') & "Backbone and SC coordinates as read from the PDB" do ires=1,nres - write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') & - ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),& + write (iout,'(i5,i3,2x,a,3f8.3,5x,3f8.3)') & + ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),& (c(j,nres+ires),j=1,3) enddo endif +! NOW LETS ROCK! SORTING + allocate(c_temporary(3,2*nres)) + allocate(itype_temporary(nres,5)) + if (.not.allocated(molnum)) allocate(molnum(nres+1)) + if (.not.allocated(istype)) write(iout,*) & + "SOMETHING WRONG WITH ISTYTPE" + allocate(istype_temp(nres)) + itype_temporary(:,:)=0 + seqalingbegin=1 + do k=1,5 + do i=1,nres + if (itype(i,k).ne.0) then + do j=1,3 + c_temporary(j,seqalingbegin)=c(j,i) + c_temporary(j,seqalingbegin+nres)=c(j,i+nres) + + enddo + itype_temporary(seqalingbegin,k)=itype(i,k) + print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin + istype_temp(seqalingbegin)=istype(i) + molnum(seqalingbegin)=k + seqalingbegin=seqalingbegin+1 + endif + enddo + enddo + do i=1,2*nres + do j=1,3 + c(j,i)=c_temporary(j,i) + enddo + enddo + do k=1,5 + do i=1,nres + itype(i,k)=itype_temporary(i,k) + istype(i)=istype_temp(i) + enddo + enddo +! if (itype(1,1).eq.ntyp1) then +! nsup=nsup-1 +! nstart_sup=2 +! if (unres_pdb) then +! 2/15/2013 by Adam: corrected insertion of the first dummy residue +! call refsys(2,3,4,e1,e2,e3,fail) +! if (fail) then +! e2(1)=0.0d0 +! e2(2)=1.0d0 +! e2(3)=0.0d0 +! endif +! do j=1,3 +! c(j,1)=c(j,2)-1.9d0*e2(j) +! enddo +! else +! do j=1,3 +! dcj=(c(j,4)-c(j,3))/2.0 +! c(j,1)=c(j,2)-dcj +! c(j,nres+1)=c(j,1) +! enddo +! endif +! endif + + if (lprn) then + write (iout,'(/a)') & + "Cartesian coordinates of the reference structure after sorting" + write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') & + "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)" + do ires=1,nres + write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') & + (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),& + (c(j,ires+nres),j=1,3) + enddo + endif +! print *,seqalingbegin,nres if(.not.allocated(vbld)) then allocate(vbld(2*nres)) do i=1,2*nres @@ -2738,11 +3556,12 @@ kkk=1 lll=0 cou=1 + write (iout,*) "symetr", symetr do i=1,nres lll=lll+1 -!c write (iout,*) "spraw lancuchy",(c(j,i),j=1,3) +! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3) if (i.gt.1) then - if ((itype(i-1).eq.ntyp1).and.(i.gt.2)) then + if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then chain_length=lll-1 kkk=kkk+1 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3) @@ -2769,7 +3588,7 @@ ! write (iout,*) "spraw lancuchy",chain_length,symetr ! do i=1,4 ! do kkk=1,chain_length -! write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3) +! write (iout,*) itype(kkk,1),(chain_rep(j,kkk,i), j=1,3) ! enddo ! enddo ! enddiagnostic @@ -2792,7 +3611,7 @@ cou=0 do kkk=1,symetr icha=tabperm(i,kkk) -! write (iout,*) i,icha + write (iout,*) i,icha do lll=1,chain_length cou=cou+1 if (cou.le.nres) then @@ -2800,7 +3619,7 @@ kupa=mod(lll,chain_length) iprzes=(kkk-1)*chain_length+lll if (kupa.eq.0) kupa=chain_length -! write (iout,*) "kupa", kupa + write (iout,*) "kupa", kupa cref(j,iprzes,i)=chain_rep(j,kupa,icha) cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha) enddo @@ -2814,16 +3633,16 @@ do kkk=1,nperm write (iout,*) "nowa struktura", nperm do i=1,nres - write (iout,110) restyp(itype(i)),i,cref(1,i,kkk),& + write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),& cref(2,i,kkk),& cref(3,i,kkk),cref(1,nres+i,kkk),& cref(2,nres+i,kkk),cref(3,nres+i,kkk) enddo - 100 format (//' alpha-carbon coordinates ',& + 100 format (//' alpha-carbon coordinates ',& ' centroid coordinates'/ & ' ', 6X,'X',11X,'Y',11X,'Z', & 10X,'X',11X,'Y',11X,'Z') - 110 format (a,'(',i3,')',6f12.5) + 110 format (a,'(',i5,')',6f12.5) enddo !c enddiag @@ -2940,9 +3759,39 @@ !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then write(iout,*) "shield_mode",shield_mode !C Varibles set size of box + with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0 + protein=index(controlcard,"PROTEIN").gt.0 + ions=index(controlcard,"IONS").gt.0 + nucleic=index(controlcard,"NUCLEIC").gt.0 + write (iout,*) "with_theta_constr ",with_theta_constr + AFMlog=(index(controlcard,'AFM')) + selfguide=(index(controlcard,'SELFGUIDE')) + print *,'AFMlog',AFMlog,selfguide,"KUPA" + call readi(controlcard,'GENCONSTR',genconstr,0) call reada(controlcard,'BOXX',boxxsize,100.0d0) call reada(controlcard,'BOXY',boxysize,100.0d0) call reada(controlcard,'BOXZ',boxzsize,100.0d0) + call readi(controlcard,'TUBEMOD',tubemode,0) + print *,"SCELE",scelemode + call readi(controlcard,"SCELEMODE",scelemode,0) + print *,"SCELE",scelemode + +! elemode = 0 is orignal UNRES electrostatics +! elemode = 1 is "Momo" potentials in progress +! elemode = 2 is in development EVALD + write (iout,*) TUBEmode,"TUBEMODE" + if (TUBEmode.gt.0) then + call reada(controlcard,"XTUBE",tubecenter(1),0.0d0) + call reada(controlcard,"YTUBE",tubecenter(2),0.0d0) + call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0) + call reada(controlcard,"RTUBE",tubeR0,0.0d0) + call reada(controlcard,"TUBETOP",bordtubetop,boxzsize) + call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0) + call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0) + buftubebot=bordtubebot+tubebufthick + buftubetop=bordtubetop-tubebufthick + endif + ! CUTOFFF ON ELECTROSTATICS call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0) call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0) @@ -3067,22 +3916,22 @@ write (iout,'(2a)') diagmeth(kdiag),& ' routine used to diagonalize matrices.' if (shield_mode.gt.0) then - pi=3.141592d0 + pi=4.0D0*datan(1.0D0) !C VSolvSphere the volume of solving sphere -!C print *,pi,"pi" + print *,pi,"pi" !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it !C there will be no distinction between proline peptide group and normal peptide !C group in case of shielding parameters - VSolvSphere=4.0/3.0*pi*rpp(1,1)**3 - VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3 + VSolvSphere=4.0/3.0*pi*(4.50d0)**3 + VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3 write (iout,*) VSolvSphere,VSolvSphere_div !C long axis of side chain - do i=1,ntyp - long_r_sidechain(i)=vbldsc0(1,i) - short_r_sidechain(i)=sigma0(i) - write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),& - sigma0(i) - enddo +! do i=1,ntyp +! long_r_sidechain(i)=vbldsc0(1,i) +! short_r_sidechain(i)=sigma0(i) +! write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),& +! sigma0(i) +! enddo buff_shield=1.0d0 endif return @@ -3207,7 +4056,7 @@ ! character(len=80) :: ucase character(len=320) :: controlcard !el local variables - integer :: i + integer :: i,j real(kind=8) :: eta call card_concat(controlcard,.true.) @@ -3242,6 +4091,10 @@ large = index(controlcard,"LARGE").gt.0 print_compon = index(controlcard,"PRINT_COMPON").gt.0 rattle = index(controlcard,"RATTLE").gt.0 + preminim=(index(controlcard,'PREMINIM').gt.0) + write (iout,*) "PREMINIM ",preminim + dccart=(index(controlcard,'CART').gt.0) + if (preminim) call read_minim ! if performing umbrella sampling, fragments constrained are read from the fragment file nset=0 if(usampl) then @@ -3320,21 +4173,27 @@ if(me.eq.king.or..not.out1file) & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",& eta - gamp=scal_fric*(pstok+rwat)*eta - stdfp=dsqrt(2*Rb*t_bath/d_time) - allocate(gamsc(ntyp1),stdfsc(ntyp1)) !(ntyp1) +! allocate(gamp + do j=1,5 !types of molecules + gamp(j)=scal_fric*(pstok(j)+rwat)*eta + stdfp(j)=dsqrt(2*Rb*t_bath/d_time) + enddo + allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1) + do j=1,5 !types of molecules do i=1,ntyp - gamsc(i)=scal_fric*(restok(i)+rwat)*eta - stdfsc(i)=dsqrt(2*Rb*t_bath/d_time) + gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta + stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time) enddo + enddo + if(me.eq.king.or..not.out1file)then write (iout,'(/2a/)') & "Radii of site types and friction coefficients and std's of",& " stochastic forces of fully exposed sites" - write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp) + write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(1),stdfp*dsqrt(gamp(1)) do i=1,ntyp - write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),& - gamsc(i),stdfsc(i)*dsqrt(gamsc(i)) + write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),& + gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1)) enddo endif else if (tbf) then @@ -3695,6 +4554,8 @@ ! Get parameter filenames and open the parameter files. call getenv_loc('BONDPAR',bondname) open (ibond,file=bondname,status='old',readonly,shared) + call getenv_loc('BONDPAR_NUCL',bondname_nucl) + open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared) call getenv_loc('THETPAR',thetname) open (ithep,file=thetname,status='old',readonly,shared) call getenv_loc('ROTPAR',rotname) @@ -3709,6 +4570,19 @@ open (ielep,file=elename,status='old',readonly,shared) call getenv_loc('SIDEPAR',sidename) open (isidep,file=sidename,status='old',readonly,shared) + + call getenv_loc('THETPAR_NUCL',thetname_nucl) + open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared) + call getenv_loc('ROTPAR_NUCL',rotname_nucl) + open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared) + call getenv_loc('TORPAR_NUCL',torname_nucl) + open (itorp_nucl,file=torname_nucl,status='old',readonly,shared) + call getenv_loc('TORDPAR_NUCL',tordname_nucl) + open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared) + call getenv_loc('SIDEPAR_NUCL',sidename_nucl) + open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared) + + #elif (defined CRAY) || (defined AIX) open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',& action='read') @@ -3719,6 +4593,9 @@ ! Get parameter filenames and open the parameter files. call getenv_loc('BONDPAR',bondname) open (ibond,file=bondname,status='old',action='read') + call getenv_loc('BONDPAR_NUCL',bondname_nucl) + open (ibond_nucl,file=bondname_nucl,status='old',action='read') + ! print *,"Processor",myrank," opened file IBOND" call getenv_loc('THETPAR',thetname) open (ithep,file=thetname,status='old',action='read') @@ -3743,8 +4620,25 @@ ! print *,"Processor",myrank," opened file IELEP" call getenv_loc('SIDEPAR',sidename) open (isidep,file=sidename,status='old',action='read') + + call getenv_loc('THETPAR_NUCL',thetname_nucl) + open (ithep_nucl,file=thetname_nucl,status='old',action='read') + call getenv_loc('ROTPAR_NUCL',rotname_nucl) + open (irotam_nucl,file=rotname_nucl,status='old',action='read') + call getenv_loc('TORPAR_NUCL',torname_nucl) + open (itorp_nucl,file=torname_nucl,status='old',action='read') + call getenv_loc('TORDPAR_NUCL',tordname_nucl) + open (itordp_nucl,file=tordname_nucl,status='old',action='read') + call getenv_loc('SIDEPAR_NUCL',sidename_nucl) + open (isidep_nucl,file=sidename_nucl,status='old',action='read') + call getenv_loc('LIPTRANPAR',liptranname) open (iliptranpar,file=liptranname,status='old',action='read') + call getenv_loc('TUBEPAR',tubename) + open (itube,file=tubename,status='old',action='read') + call getenv_loc('IONPAR',ionname) + open (iion,file=ionname,status='old',action='read') + ! print *,"Processor",myrank," opened file ISIDEP" ! print *,"Processor",myrank," opened parameter files" #elif (defined G77) @@ -3754,6 +4648,9 @@ ! Get parameter filenames and open the parameter files. call getenv_loc('BONDPAR',bondname) open (ibond,file=bondname,status='old') + call getenv_loc('BONDPAR_NUCL',bondname_nucl) + open (ibond_nucl,file=bondname_nucl,status='old') + call getenv_loc('THETPAR',thetname) open (ithep,file=thetname,status='old') call getenv_loc('ROTPAR',rotname) @@ -3770,8 +4667,23 @@ open (ielep,file=elename,status='old') call getenv_loc('SIDEPAR',sidename) open (isidep,file=sidename,status='old') + + open (ithep_nucl,file=thetname_nucl,status='old') + call getenv_loc('ROTPAR_NUCL',rotname_nucl) + open (irotam_nucl,file=rotname_nucl,status='old') + call getenv_loc('TORPAR_NUCL',torname_nucl) + open (itorp_nucl,file=torname_nucl,status='old') + call getenv_loc('TORDPAR_NUCL',tordname_nucl) + open (itordp_nucl,file=tordname_nucl,status='old') + call getenv_loc('SIDEPAR_NUCL',sidename_nucl) + open (isidep_nucl,file=sidename_nucl,status='old') + call getenv_loc('LIPTRANPAR',liptranname) open (iliptranpar,file=liptranname,status='old') + call getenv_loc('TUBEPAR',tubename) + open (itube,file=tubename,status='old') + call getenv_loc('IONPAR',ionname) + open (iion,file=ionname,status='old') #else open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',& readonly) @@ -3780,6 +4692,8 @@ ! Get parameter filenames and open the parameter files. call getenv_loc('BONDPAR',bondname) open (ibond,file=bondname,status='old',action='read') + call getenv_loc('BONDPAR_NUCL',bondname_nucl) + open (ibond_nucl,file=bondname_nucl,status='old',action='read') call getenv_loc('THETPAR',thetname) open (ithep,file=thetname,status='old',action='read') call getenv_loc('ROTPAR',rotname) @@ -3802,13 +4716,50 @@ open (ielep,file=elename,status='old',readonly) call getenv_loc('SIDEPAR',sidename) open (isidep,file=sidename,status='old',readonly) + + call getenv_loc('THETPAR_NUCL',thetname_nucl) + open (ithep_nucl,file=thetname_nucl,status='old',action='read') + call getenv_loc('ROTPAR_NUCL',rotname_nucl) + open (irotam_nucl,file=rotname_nucl,status='old',action='read') + call getenv_loc('TORPAR_NUCL',torname_nucl) + open (itorp_nucl,file=torname_nucl,status='old',action='read') + call getenv_loc('TORDPAR_NUCL',tordname_nucl) + open (itordp_nucl,file=tordname_nucl,status='old',action='read') + call getenv_loc('SIDEPAR_NUCL',sidename_nucl) + open (isidep_nucl,file=sidename_nucl,status='old',action='read') + call getenv_loc('SIDEPAR_SCBASE',sidename_scbase) + open (isidep_scbase,file=sidename_scbase,status='old',action='read') + call getenv_loc('PEPPAR_PEPBASE',pepname_pepbase) + open (isidep_pepbase,file=pepname_pepbase,status='old',action='read') + call getenv_loc('SCPAR_PHOSPH',pepname_scpho) + open (isidep_scpho,file=pepname_scpho,status='old',action='read') + call getenv_loc('PEPPAR_PHOSPH',pepname_peppho) + open (isidep_peppho,file=pepname_peppho,status='old',action='read') + + call getenv_loc('LIPTRANPAR',liptranname) open (iliptranpar,file=liptranname,status='old',action='read') + call getenv_loc('TUBEPAR',tubename) + open (itube,file=tubename,status='old',action='read') + call getenv_loc('IONPAR',ionname) + open (iion,file=ionname,status='old',action='read') + #ifndef CRYST_SC call getenv_loc('ROTPARPDB',rotname_pdb) open (irotam_pdb,file=rotname_pdb,status='old',action='read') #endif #endif + call getenv_loc('SCPPAR_NUCL',scpname_nucl) +#if defined(WINIFL) || defined(WINPGI) + open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared) +#elif (defined CRAY) || (defined AIX) + open (iscpp_nucl,file=scpname_nucl,status='old',action='read') +#elif (defined G77) + open (iscpp_nucl,file=scpname_nucl,status='old') +#else + open (iscpp_nucl,file=scpname_nucl,status='old',action='read') +#endif + #ifndef OLDSCP ! ! 8/9/01 In the newest version SCp interaction constants are read from a file @@ -4000,6 +4951,7 @@ open(irest2,file=rest2name,status='unknown') read(irest2,*) totT,EK,potE,totE,t_bath + totTafm=totT ! do i=1,2*nres ! AL 4/17/17: Now reading d_t(0,:) too do i=0,2*nres