! 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 ',&
! 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) ='-'
ii=0
do i=1,nres
+ ftors(i)=ftors(1)
if ( secstruc(i) .eq. 'H') then
! Helix restraints for this residue
ii=ii+1
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))
!
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))
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
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
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)
' 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-))') &
' 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-))') &
' 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
' 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
' 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-))') &
' 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
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)
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)
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.
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
#endif
if (lprint) then
+ l=3
write (iout,'(/a/)') 'Torsional constants:'
do i=1,nsccortyp
do j=1,nsccortyp
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
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))
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 --------------------------------
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
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)
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
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
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)
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 ',&
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
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)
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)
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
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
!
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
! 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'
! 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
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.
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
! 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
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
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
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
! 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
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
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
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
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
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)
! 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
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
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
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
!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)
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
! 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.)
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
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
! 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)
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')
! 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')
! 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)
! 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)
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)
! 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)
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
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
projects / unres4.git / blobdiff