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,&
!
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))
endif
enddo
#else
- read (ibond,*) junk,vbldp0,vbldpDUM,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,1),nbondterm(i),&
- vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(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)
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)
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))
end select
continue
close (isidep)
+
!-----------------------------------------------------------------------
! Calculate the "working" parameters of SC interactions.
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
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)))
+
+! 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)
+
+ 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
!
! 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
character(len=80) :: card
real(kind=8),dimension(3,20) :: sccor
integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode,
- integer :: isugar
+ integer :: isugar,molecprev,firstion
character*1 :: sugar
real(kind=8) :: cou
real(kind=8),dimension(3) :: ccc
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
ishift1=ishift1-1 !!!!!
! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
ires=ires-ishift+ishift1
- print *,ires,ishift,ishift1
+! print *,ires,ishift,ishift1
ires_old=ires
ibeg=0
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
molecule=2
- itype(ires,molecule)=rescode(ires,res(3:4),0,molecule)
+ 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
.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
+! print *,ires,ishift,ishift1
counter=counter+1
! iii=iii+1
! do j=1,3
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(1:3).eq."C5'") then
+ 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)
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 ') then
+ 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
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) nres_molec(molecule)=nres_molec(molecule)-2
+ if ((ires_old.ne.ires).and.(molecule.ne.5)) &
+ nres_molec(molecule)=nres_molec(molecule)-2
! print *,'I have', nres_molec(:)
- do k=1,5
+ 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,1)
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"
+! print *,"masz sieczke"
do j=1,5
if (itype(i+2,j).ne.0) then
itype(i+1,k)=0
! 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
+! 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
+! endif !unres_pdb
else !itype(i+1,1).eq.ntyp1
- if (unres_pdb) then
+! 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 !unres_pdb
endif !itype(i+1,1).eq.ntyp1
endif !itype.eq.ntyp1
nsup=nres
nstart_sup=1
! print *,"molecule",molecule
- if (itype(nres,1).ne.10) then
+ if ((itype(nres,1).ne.10)) then
nres=nres+1
+ if (molecule.eq.5) molecule=molecprev
itype(nres,molecule)=ntyp1_molec(molecule)
nres_molec(molecule)=nres_molec(molecule)+1
- if (unres_pdb) then
+! 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(:)
! NOW LETS ROCK! SORTING
allocate(c_temporary(3,2*nres))
allocate(itype_temporary(nres,5))
+ allocate(molnum(nres+1))
+ allocate(istype_temp(nres))
itype_temporary(:,:)=0
seqalingbegin=1
do k=1,5
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
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
enddo
endif
- print *,seqalingbegin,nres
+! 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)
! 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,1),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"
! 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
projects / unres4.git / blobdiff