! 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 ',&
use geometry_data
use energy_data
- use control_data, only:maxtor,maxterm
+ use control_data, only:maxterm !,maxtor
use MD_data
use MPI_data
!el use map_data
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
+ 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(msc(ntyp+1,5)) !(ntyp+1)
+ allocate(isc(ntyp+1,5)) !(ntyp+1)
+ allocate(restok(ntyp+1,5)) !(ntyp+1)
allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
-
+ allocate(long_r_sidechain(ntyp))
+ allocate(short_r_sidechain(ntyp))
dsc(:)=0.0d0
dsc_inv(:)=0.0d0
endif
enddo
#else
- read (ibond,*) junk,vbldp0,akp,rjunk,mp,ip,pstok
- do i=1,ntyp
+ 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)
allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
- do i=-ntyp,ntyp
- a0thet(i)=0.0D0
- do j=1,2
- do ichir1=-1,1
- do ichir2=-1,1
- athet(j,i,ichir1,ichir2)=0.0D0
- bthet(j,i,ichir1,ichir2)=0.0D0
- enddo
- enddo
- enddo
- do j=0,3
- polthet(j,i)=0.0D0
- enddo
- do j=1,3
- gthet(j,i)=0.0D0
- enddo
- theta0(i)=0.0D0
- sig0(i)=0.0D0
- sigc0(i)=0.0D0
- enddo
+
+ a0thet(:)=0.0D0
+ athet(:,:,:,:)=0.0D0
+ bthet(:,:,:,:)=0.0D0
+ polthet(:,:)=0.0D0
+ gthet(:,:)=0.0D0
+ theta0(:)=0.0D0
+ sig0(:)=0.0D0
+ sigc0(:)=0.0D0
+ allocate(liptranene(ntyp))
+!C reading lipid parameters
+ write (iout,*) "iliptranpar",iliptranpar
+ call flush(iout)
+ read(iliptranpar,*) pepliptran
+ print *,pepliptran
+ do i=1,ntyp
+ read(iliptranpar,*) liptranene(i)
+ print *,liptranene(i)
+ enddo
+ close(iliptranpar)
#ifdef CRYST_THETA
!
' 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
!----------------------------------------------------
allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
- allocate(aa0thet(-maxthetyp1:maxthetyp1,&
- -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
+ allocate(aa0thet(-nthetyp-1:nthetyp+1,&
+ -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
!(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
- allocate(aathet(ntheterm,-maxthetyp1:maxthetyp1,&
- -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
+ allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
+ -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
!(maxtheterm,-maxthetyp1:maxthetyp1,&
! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
- allocate(bbthet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,&
- -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
- allocate(ccthet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,&
- -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
- allocate(ddthet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,&
- -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
- allocate(eethet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,&
- -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
+ allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
+ -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
+ allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
+ -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
+ allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
+ -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
+ allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
+ -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
!(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
- allocate(ffthet(ndouble,ndouble,ntheterm3,-maxthetyp1:maxthetyp1,&
- -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
- allocate(ggthet(ndouble,ndouble,ntheterm3,-maxthetyp1:maxthetyp1,&
- -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
+ allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
+ -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
+ allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
+ -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
!(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
do i=-ntyp1,-1
ithetyp(i)=-ithetyp(-i)
enddo
- do iblock=1,2
- do i=-maxthetyp,maxthetyp
- do j=-maxthetyp,maxthetyp
- do k=-maxthetyp,maxthetyp
- aa0thet(i,j,k,iblock)=0.0d0
- do l=1,ntheterm
- aathet(l,i,j,k,iblock)=0.0d0
- enddo
- do l=1,ntheterm2
- do m=1,nsingle
- bbthet(m,l,i,j,k,iblock)=0.0d0
- ccthet(m,l,i,j,k,iblock)=0.0d0
- ddthet(m,l,i,j,k,iblock)=0.0d0
- eethet(m,l,i,j,k,iblock)=0.0d0
- enddo
- enddo
- do l=1,ntheterm3
- do m=1,ndouble
- do mm=1,ndouble
- ffthet(mm,m,l,i,j,k,iblock)=0.0d0
- ggthet(mm,m,l,i,j,k,iblock)=0.0d0
- enddo
- enddo
- enddo
- enddo
- enddo
- enddo
- enddo
+
+ aa0thet(:,:,:,:)=0.0d0
+ aathet(:,:,:,:,:)=0.0d0
+ bbthet(:,:,:,:,:,:)=0.0d0
+ ccthet(:,:,:,:,:,:)=0.0d0
+ ddthet(:,:,:,:,:,:)=0.0d0
+ eethet(:,:,:,:,:,:)=0.0d0
+ ffthet(:,:,:,:,:,:,:)=0.0d0
+ ggthet(:,:,:,:,:,:,:)=0.0d0
+
! VAR:iblock means terminally blocking group 1=non-proline 2=proline
do iblock=1,2
! VAR:ntethtyp is type of theta potentials type currently 0=glycine
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)
allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
- do i=1,ntyp
- do j=1,maxlob
- bsc(j,i)=0.0D0
- nlob(i)=0
- enddo
- enddo
- nlob(ntyp1)=0
- dsc(ntyp1)=0.0D0
-
- do i=-ntyp,ntyp
- do j=1,maxlob
- do k=1,3
- censc(k,j,i)=0.0D0
- enddo
- do k=1,3
- do l=1,3
- gaussc(l,k,j,i)=0.0D0
- enddo
- enddo
- enddo
- enddo
-
+ bsc(:,:)=0.0D0
+ nlob(:)=0
+ nlob(:)=0
+ dsc(:)=0.0D0
+ censc(:,:,:)=0.0D0
+ gaussc(:,:,:,:)=0.0D0
+
#ifdef CRYST_SC
!
! Read the parameters of the probability distribution/energy expression
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.
allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
!el---------------------------
- do iblock=1,2
- do i=-ntortyp,ntortyp
- do j=-ntortyp,ntortyp
- nterm(i,j,iblock)=0
- nlor(i,j,iblock)=0
- enddo
- enddo
- enddo
+ nterm(:,:,:)=0
+ nlor(:,:,:)=0
!el---------------------------
read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
do i=-ntyp,-1
itortyp(i)=-itortyp(-i)
enddo
-! itortyp(ntyp1)=ntortyp
-! itortyp(-ntyp1)=-ntortyp
+ itortyp(ntyp1)=ntortyp
+ itortyp(-ntyp1)=-ntortyp
do iblock=1,2
write (iout,*) 'ntortyp',ntortyp
do i=0,ntortyp-1
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
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)
- do i=1,ntyp
- do j=1,ntyp
- augm(i,j)=0.0D0
- enddo
- chip(i)=0.0D0
- alp(i)=0.0D0
- sigma0(i)=0.0D0
- sigii(i)=0.0D0
- rr0(i)=0.0D0
- enddo
+ allocate(epslip(ntyp,ntyp))
+ augm(:,:)=0.0D0
+ chip(:)=0.0D0
+ alp(:)=0.0D0
+ sigma0(:)=0.0D0
+ sigii(:)=0.0D0
+ rr0(:)=0.0D0
+
!--------------------------------
read (isidep,*,end=117,err=117) ipot,expon
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
read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
+ do i=1,ntyp
+ read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
+ enddo
+
! For the GB potential convert sigma'**2 into chi'
if (ipot.eq.4) then
do i=1,ntyp
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
! goto 50
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(ntyp1,ntyp1),bb(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
+ allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1),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)
- do i=1,ntyp1
- do j=1,ntyp1
- aa(i,j)=0.0D0
- bb(i,j)=0.0D0
- chi(i,j)=0.0D0
- sigma(i,j)=0.0D0
- r0(i,j)=0.0D0
- enddo
- enddo
+ 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
+ chi(:,:)=0.0D0
+ sigma(:,:)=0.0D0
+ r0(:,:)=0.0D0
+ 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
do j=1,i-1
eps(i,j)=eps(j,i)
+ epslip(i,j)=epslip(j,i)
enddo
enddo
do i=1,ntyp
epsij=eps(i,j)
sigeps=dsign(1.0D0,epsij)
epsij=dabs(epsij)
- aa(i,j)=epsij*rrij*rrij
- bb(i,j)=-sigeps*epsij*rrij
- aa(j,i)=aa(i,j)
- bb(j,i)=bb(i,j)
+ aa_aq(i,j)=epsij*rrij*rrij
+ bb_aq(i,j)=-sigeps*epsij*rrij
+ aa_aq(j,i)=aa_aq(i,j)
+ bb_aq(j,i)=bb_aq(i,j)
+ epsijlip=epslip(i,j)
+ sigeps=dsign(1.0D0,epsijlip)
+ epsijlip=dabs(epsijlip)
+ aa_lip(i,j)=epsijlip*rrij*rrij
+ bb_lip(i,j)=-sigeps*epsijlip*rrij
+ 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
sigt1sq=sigma0(i)**2
sigt2sq=sigma0(j)**2
endif
if (lprint) then
write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
- restyp(i),restyp(j),aa(i,j),bb(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(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
- do i=1,ntyp
- do j=1,2
- bad(i,j)=0.0D0
+ 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
+
+ allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
+ bad(:,:)=0.0D0
+
#ifdef OLDSCP
!
! Define the SC-p interaction constants (hard-coded; old style)
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
!
! 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
goto 10
else if (card(:3).eq.'TER') then
! End current chain
- ires_old=ires+1
+ ires_old=ires+2
+ ishift=ishift+1
ishift1=ishift1+1
- itype(ires_old)=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
enddo
else
call sccenter(ires,iii,sccor)
+! iii=0
endif
iii=0
endif
! 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
! write (iout,*) "Calculating sidechain center iii",iii
if (unres_pdb) then
do j=1,3
- dc(j,ires+nres)=sccor(j,iii)
+ dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
enddo
else
call sccenter(ires_old,iii,sccor)
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
- itype(ires)=rescode(ires,res,0)
+ 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
+ molecule=2
+ itype(ires,molecule)=rescode(ires,res(2:4),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
+ read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
+ else
+ iii=iii+1
+ read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
+ endif
+! 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)
+ itype(ires,molecule)=rescode(ires,res(2:4),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_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)
- 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
+! 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
+! dc(j,i)=c(j,i)
+! enddo
+! endif
+ 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,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
+ 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,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
+ 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+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,1).eq.ntyp1
+ endif !itype.eq.ntyp1
+
+ enddo
enddo
! Calculate the CM of the last side chain.
if (iii.gt.0) 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 (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)
e2(3)=0.0d0
endif
do j=1,3
- c(j,nres)=c(j,nres-1)-3.8d0*e2(j)
+ 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)
+ 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
+! 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
+ write (iout,*) "Error: wrong parameter value: NRES=",nres,&
+ " NRES0=",nres0
+ stop "Error nres value in WHAM input"
+ endif
+#endif
!---------------------------------
!el reallocate tables
! do i=1,maxres/3
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
e2(3)=0.0d0
endif
do j=1,3
- c(j,1)=c(j,2)-3.8d0*e2(j)
+ c(j,1)=c(j,2)-1.9d0*e2(j)
enddo
else
do j=1,3
- dcj=c(j,4)-c(j,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
+! 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
write (iout,'(a,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),i3,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
+! znamy już nres wiec mozna alokowac tablice
! Calculate internal coordinates.
if(me.eq.king.or..not.out1file)then
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),&
+ 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))
+ allocate(molnum(nres))
+ 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,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),i3,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
!!!el
if(.not.allocated(theta)) then
allocate(theta(nres+2))
-! allocate(phi(nres+2))
-! allocate(alph(nres+2))
-! allocate(omeg(nres+2))
- do i=1,nres+2
- theta(i)=0.0d0
-! phi(i)=0.0d0
-! alph(i)=0.0d0
-! omeg(i)=0.0d0
- enddo
+ theta(:)=0.0d0
endif
-! allocate(costtab(nres))
-! allocate(sinttab(nres))
-! allocate(cost2tab(nres))
-! allocate(sint2tab(nres))
-! allocate(xxref(nres))
-! allocate(yyref(nres))
-! allocate(zzref(nres)) !(maxres)
-! do i=1,nres
-! costtab(i)=0.0d0
-! sinttab(i)=0.0d0
-! cost2tab(i)=0.0d0
-! sint2tab(i)=0.0d0
-! xxref(i)=0.0d0
-! yyref(i)=0.0d0
-! zzref(i)=0.0d0
-! enddo
-
-! endif
+
if(.not.allocated(phi)) allocate(phi(nres+2))
if(.not.allocated(alph)) allocate(alph(nres+2))
if(.not.allocated(omeg)) allocate(omeg(nres+2))
if(.not.allocated(dc_norm)) then
! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
allocate(dc_norm(3,0:2*nres+2))
- do i=0,2*nres+2
- dc_norm(1,i)=0.d0
- dc_norm(2,i)=0.d0
- dc_norm(3,i)=0.d0
- enddo
+ dc_norm(:,:)=0.d0
endif
call int_from_cart(.true.,.false.)
- call sc_loc_geom(.true.)
-! call sc_loc_geom(.false.)
-! wczesbiej bylo false
+ call sc_loc_geom(.false.)
do i=1,nres
thetaref(i)=theta(i)
phiref(i)=phi(i)
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)
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)
! enddiagnostic
! makes copy of chains
write (iout,*) "symetr", symetr
+ do j=1,3
+ dc(j,0)=c(j,1)
+ enddo
if (symetr.gt.1) then
call permut(symetr)
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)
character(len=640) :: controlcard
real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
-
+ integer i
nglob_csa=0
eglob_csa=1d99
timem=timlim
modecalc=0
call reada(controlcard,"T_BATH",t_bath,300.0d0)
+!C SHIELD keyword sets if the shielding effect of side-chains is used
+!C 0 denots no shielding is used all peptide are equally despite the
+!C solvent accesible area
+!C 1 the newly introduced function
+!C 2 reseved for further possible developement
+ call readi(controlcard,'SHIELD',shield_mode,0)
+!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)
+ 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,*) "R_CUT_ELE=",r_cut_ele
+! Lipidic parameters
+ call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
+ call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
+ if (lipthick.gt.0.0d0) then
+ bordliptop=(boxzsize+lipthick)/2.0
+ bordlipbot=bordliptop-lipthick
+ if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
+ write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
+ buflipbot=bordlipbot+lipbufthick
+ bufliptop=bordliptop-lipbufthick
+ if ((lipbufthick*2.0d0).gt.lipthick) &
+ write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
+ endif !lipthick.gt.0
+ write(iout,*) "bordliptop=",bordliptop
+ write(iout,*) "bordlipbot=",bordlipbot
+ write(iout,*) "bufliptop=",bufliptop
+ write(iout,*) "buflipbot=",buflipbot
+ write (iout,*) "SHIELD MODE",shield_mode
+
+!C-------------------------
minim=(index(controlcard,'MINIMIZE').gt.0)
dccart=(index(controlcard,'CART').gt.0)
overlapsc=(index(controlcard,'OVERLAP').gt.0)
if(me.eq.king.or..not.out1file) &
write (iout,'(2a)') diagmeth(kdiag),&
' routine used to diagonalize matrices.'
+ if (shield_mode.gt.0) then
+ pi=3.141592d0
+!C VSolvSphere the volume of solving sphere
+!C 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
+ 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
+ buff_shield=1.0d0
+ endif
return
end subroutine read_control
!-----------------------------------------------------------------------------
! 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.)
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
!-----------------------------------------------------------------------------
subroutine openunits
- use energy_data, only: usampl
+ use MD_data, only: usampl
use csa_data
use MPI_data
use control_data, only:out1file
! 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('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
subroutine readrst
use geometry_data, only: nres,dc
- use energy_data, only: usampl,iset
use MD_data
! implicit real*8 (a-h,o-z)
! include 'DIMENSIONS'
open(irest2,file=rest2name,status='unknown')
read(irest2,*) totT,EK,potE,totE,t_bath
- do i=1,2*nres
+ totTafm=totT
+! do i=1,2*nres
+! AL 4/17/17: Now reading d_t(0,:) too
+ do i=0,2*nres
read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
enddo
- do i=1,2*nres
+! do i=1,2*nres
+! AL 4/17/17: Now reading d_c(0,:) too
+ do i=0,2*nres
read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
enddo
if(usampl) then