include 'COMMON.SBRIDGE'
include 'COMMON.MD'
include 'COMMON.SETUP'
+ include 'COMMON.CONTROL'
+ include 'COMMON.SHIELD'
character*1 t1,t2,t3
character*1 onelett(4) /"G","A","P","D"/
+ character*1 toronelet(-2:2) /"p","a","G","A","P"/
logical lprint,LaTeX
dimension blower(3,3,maxlob)
- dimension b(13)
+ character*3 string
+C dimension b(13)
character*3 lancuch,ucase
C
C For printing parameters after they are read set the following in the UNRES
vblinv2=vblinv*vblinv
c
c Read the virtual-bond parameters, masses, and moments of inertia
-c and Stokes' radii of the peptide group and side chains
+c and Stokes radii of the peptide group and side chains
c
#ifdef CRYST_BOND
- read (ibond,*) vbldp0,akp,mp,ip,pstok
+ read (ibond,*) vbldp0,vbldpdum,akp,mp,ip,pstok
do i=1,ntyp
nbondterm(i)=1
read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
endif
enddo
#else
- read (ibond,*) junk,vbldp0,akp,rjunk,mp,ip,pstok
+ read (ibond,*) junk,vbldp0,vbldpdum,akp,rjunk,mp,ip,pstok
do i=1,ntyp
read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
& j=1,nbondterm(i)),msc(i),isc(i),restok(i)
enddo
enddo
endif
+C reading lipid parameters
+ if (lprint) then
+ write (iout,*) "iliptranpar",iliptranpar
+ call flush(iout)
+ endif
+ read(iliptranpar,*) pepliptran
+ do i=1,ntyp
+ read(iliptranpar,*) liptranene(i)
+ enddo
+ close(iliptranpar)
#ifdef CRYST_THETA
C
C Read the parameters of the probability distribution/energy expression
C of the virtual-bond valence angles theta
C
do i=1,ntyp
- read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i),j=1,2),
- & (bthet(j,i),j=1,2)
+ read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),
+ & (bthet(j,i,1,1),j=1,2)
read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
- read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
- read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
- sigc0(i)=sigc0(i)**2
+ read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
+ read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
+ sigc0(i)=sigc0(i)**2
+ enddo
+ do i=1,ntyp
+ athet(1,i,1,-1)=athet(1,i,1,1)
+ athet(2,i,1,-1)=athet(2,i,1,1)
+ bthet(1,i,1,-1)=-bthet(1,i,1,1)
+ bthet(2,i,1,-1)=-bthet(2,i,1,1)
+ athet(1,i,-1,1)=-athet(1,i,1,1)
+ athet(2,i,-1,1)=-athet(2,i,1,1)
+ bthet(1,i,-1,1)=bthet(1,i,1,1)
+ bthet(2,i,-1,1)=bthet(2,i,1,1)
+ enddo
+ do i=-ntyp,-1
+ a0thet(i)=a0thet(-i)
+ athet(1,i,-1,-1)=athet(1,-i,1,1)
+ athet(2,i,-1,-1)=-athet(2,-i,1,1)
+ bthet(1,i,-1,-1)=bthet(1,-i,1,1)
+ bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
+ athet(1,i,-1,1)=athet(1,-i,1,1)
+ athet(2,i,-1,1)=-athet(2,-i,1,1)
+ bthet(1,i,-1,1)=-bthet(1,-i,1,1)
+ bthet(2,i,-1,1)=bthet(2,-i,1,1)
+ athet(1,i,1,-1)=-athet(1,-i,1,1)
+ athet(2,i,1,-1)=athet(2,-i,1,1)
+ bthet(1,i,1,-1)=bthet(1,-i,1,1)
+ bthet(2,i,1,-1)=-bthet(2,-i,1,1)
+ theta0(i)=theta0(-i)
+ sig0(i)=sig0(-i)
+ sigc0(i)=sigc0(-i)
+ do j=0,3
+ polthet(j,i)=polthet(j,-i)
+ enddo
+ do j=1,3
+ gthet(j,i)=gthet(j,-i)
+ enddo
enddo
+
close (ithep)
if (lprint) then
if (.not.LaTeX) then
& ' B1 ',' B2 '
do i=1,ntyp
write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
- & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
+ & 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-))')
& 'Parameters of the expression for sigma(theta_c):',
& ' b1*10^1 ',' b2*10^1 '
do i=1,ntyp
write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
- & a0thet(i),(100*athet(j,i),j=1,2),(10*bthet(j,i),j=1,2)
+ & a0thet(i),(100*athet(j,i,1,1),j=1,2),
+ & (10*bthet(j,i,1,1),j=1,2)
enddo
write (iout,'(/a/9x,5a/79(1h-))')
& 'Parameters of the expression for sigma(theta_c):',
endif
endif
#else
+ IF (tor_mode.eq.0) THEN
C
C Read the parameters of Utheta determined from ab initio surfaces
C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
C
read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,
& ntheterm3,nsingle,ndouble
+ write (iout,*) "ithep",ithep
+ call flush(iout)
nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
- do i=1,maxthetyp
- do j=1,maxthetyp
- do k=1,maxthetyp
- aa0thet(i,j,k)=0.0d0
+ 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)=0.0d0
+ aathet(l,i,j,k,iblock)=0.0d0
enddo
do l=1,ntheterm2
do m=1,nsingle
- bbthet(m,l,i,j,k)=0.0d0
- ccthet(m,l,i,j,k)=0.0d0
- ddthet(m,l,i,j,k)=0.0d0
- eethet(m,l,i,j,k)=0.0d0
+ 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)=0.0d0
- ggthet(mm,m,l,i,j,k)=0.0d0
+ 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
- do i=1,nthetyp
- do j=1,nthetyp
- do k=1,nthetyp
- read (ithep,'(3a)',end=111,err=111) res1,res2,res3
- read (ithep,*,end=111,err=111) aa0thet(i,j,k)
- read (ithep,*,end=111,err=111)(aathet(l,i,j,k),l=1,ntheterm)
+ enddo
+ enddo
+c VAR:iblock means terminally blocking group 1=non-proline 2=proline
+ do iblock=1,2
+c VAR:ntethtyp is type of theta potentials type currently 0=glycine
+c VAR:1=non-glicyne non-proline 2=proline
+c VAR:negative values for D-aminoacid
+ do i=0,nthetyp
+ do j=-nthetyp,nthetyp
+ do k=-nthetyp,nthetyp
+ read (ithep,'(6a)',end=111,err=111) res1
+ read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
+c VAR: aa0thet is variable describing the average value of Foureir
+c VAR: expansion series
+c VAR: aathet is foureir expansion in theta/2 angle for full formula
+c VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
+Condens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
+ read (ithep,*,end=111,err=111)
+ &(aathet(l,i,j,k,iblock),l=1,ntheterm)
read (ithep,*,end=111,err=111)
- & ((bbthet(lll,ll,i,j,k),lll=1,nsingle),
- & (ccthet(lll,ll,i,j,k),lll=1,nsingle),
- & (ddthet(lll,ll,i,j,k),lll=1,nsingle),
- & (eethet(lll,ll,i,j,k),lll=1,nsingle),ll=1,ntheterm2)
+ & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
+ & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
+ & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
+ & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),
+ & ll=1,ntheterm2)
read (ithep,*,end=111,err=111)
- & (((ffthet(llll,lll,ll,i,j,k),ffthet(lll,llll,ll,i,j,k),
- & ggthet(llll,lll,ll,i,j,k),ggthet(lll,llll,ll,i,j,k),
+ & (((ffthet(llll,lll,ll,i,j,k,iblock),
+ & ffthet(lll,llll,ll,i,j,k,iblock),
+ & ggthet(llll,lll,ll,i,j,k,iblock),
+ & ggthet(lll,llll,ll,i,j,k,iblock),
& llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
enddo
enddo
C
C For dummy ends assign glycine-type coefficients of theta-only terms; the
C coefficients of theta-and-gamma-dependent terms are zero.
-C
+C IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
+C RECOMENTDED AFTER VERSION 3.3)
+c do i=1,nthetyp
+c do j=1,nthetyp
+c do l=1,ntheterm
+c aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
+c aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
+c enddo
+c aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
+c aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
+c enddo
+c do l=1,ntheterm
+c aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
+c enddo
+c aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
+c enddo
+c enddo
+C AND COMMENT THE LOOPS BELOW
do i=1,nthetyp
do j=1,nthetyp
do l=1,ntheterm
- aathet(l,i,j,nthetyp+1)=aathet(l,i,j,1)
- aathet(l,nthetyp+1,i,j)=aathet(l,1,i,j)
+ aathet(l,i,j,nthetyp+1,iblock)=0.0d0
+ aathet(l,nthetyp+1,i,j,iblock)=0.0d0
enddo
- aa0thet(i,j,nthetyp+1)=aa0thet(i,j,1)
- aa0thet(nthetyp+1,i,j)=aa0thet(1,i,j)
+ aa0thet(i,j,nthetyp+1,iblock)=0.0d0
+ aa0thet(nthetyp+1,i,j,iblock)=0.0d0
enddo
do l=1,ntheterm
- aathet(l,nthetyp+1,i,nthetyp+1)=aathet(l,1,i,1)
+ aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
enddo
- aa0thet(nthetyp+1,i,nthetyp+1)=aa0thet(1,i,1)
+ aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
+ enddo
enddo
+C TILL HERE
+C Substitution for D aminoacids from symmetry.
+ do iblock=1,2
+ do i=-nthetyp,0
+ do j=-nthetyp,nthetyp
+ do k=-nthetyp,nthetyp
+ aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
+ do l=1,ntheterm
+ aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
+ enddo
+ do ll=1,ntheterm2
+ do lll=1,nsingle
+ bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
+ ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
+ ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
+ eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
+ enddo
+ enddo
+ do ll=1,ntheterm3
+ do lll=2,ndouble
+ do llll=1,lll-1
+ ffthet(llll,lll,ll,i,j,k,iblock)=
+ & ffthet(llll,lll,ll,-i,-j,-k,iblock)
+ ffthet(lll,llll,ll,i,j,k,iblock)=
+ & ffthet(lll,llll,ll,-i,-j,-k,iblock)
+ ggthet(llll,lll,ll,i,j,k,iblock)=
+ & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
+ ggthet(lll,llll,ll,i,j,k,iblock)=
+ & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
+ enddo !ll
+ enddo !lll
+ enddo !llll
+ enddo !k
+ enddo !j
+ enddo !i
+ enddo !iblock
C
C Control printout of the coefficients of virtual-bond-angle potentials
C
if (lprint) then
write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
- do i=1,nthetyp+1
- do j=1,nthetyp+1
- do k=1,nthetyp+1
+ do iblock=1,2
+ do i=0,nthetyp
+ do j=-nthetyp,nthetyp
+ do k=-nthetyp,nthetyp
write (iout,'(//4a)')
- & 'Type ',onelett(i),onelett(j),onelett(k)
+ & 'Type ',toronelet(i),toronelet(j),toronelet(k)
write (iout,'(//a,10x,a)') " l","a[l]"
- write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k)
+ write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
write (iout,'(i2,1pe15.5)')
- & (l,aathet(l,i,j,k),l=1,ntheterm)
+ & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
do l=1,ntheterm2
write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))')
& "b",l,"c",l,"d",l,"e",l
do m=1,nsingle
write (iout,'(i2,4(1pe15.5))') m,
- & bbthet(m,l,i,j,k),ccthet(m,l,i,j,k),
- & ddthet(m,l,i,j,k),eethet(m,l,i,j,k)
+ & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
+ & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
enddo
enddo
do l=1,ntheterm3
do m=2,ndouble
do n=1,m-1
write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
- & ffthet(n,m,l,i,j,k),ffthet(m,n,l,i,j,k),
- & ggthet(n,m,l,i,j,k),ggthet(m,n,l,i,j,k)
+ & ffthet(n,m,l,i,j,k,iblock),
+ & ffthet(m,n,l,i,j,k,iblock),
+ & ggthet(n,m,l,i,j,k,iblock),
+ & ggthet(m,n,l,i,j,k,iblock)
enddo
enddo
enddo
enddo
enddo
+ enddo
enddo
call flush(iout)
endif
+
+ ELSE
+
+C here will be the apropriate recalibrating for D-aminoacid
+ read (ithep,*,end=121,err=121) nthetyp
+ do i=-nthetyp+1,nthetyp-1
+ read (ithep,*,end=121,err=121) nbend_kcc_Tb(i)
+ do j=0,nbend_kcc_Tb(i)
+ read (ithep,*,end=121,err=121) ijunk,v1bend_chyb(j,i)
+ enddo
+ enddo
+ if (lprint) then
+ write (iout,'(a)')
+ & "Parameters of the valence-only potentials"
+ do i=-nthetyp+1,nthetyp-1
+ write (iout,'(2a)') "Type ",toronelet(i)
+ do k=0,nbend_kcc_Tb(i)
+ write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
+ enddo
+ enddo
+ endif
+
+ ENDIF ! TOR_MODE
+
+c write (2,*) "Start reading THETA_PDB",ithep_pdb
+ do i=1,ntyp
+c write (2,*) 'i=',i
+ read (ithep_pdb,*,err=111,end=111)
+ & a0thet(i),(athet(j,i,1,1),j=1,2),
+ & (bthet(j,i,1,1),j=1,2)
+ read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
+ read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
+ read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
+ sigc0(i)=sigc0(i)**2
+ enddo
+ do i=1,ntyp
+ athet(1,i,1,-1)=athet(1,i,1,1)
+ athet(2,i,1,-1)=athet(2,i,1,1)
+ bthet(1,i,1,-1)=-bthet(1,i,1,1)
+ bthet(2,i,1,-1)=-bthet(2,i,1,1)
+ athet(1,i,-1,1)=-athet(1,i,1,1)
+ athet(2,i,-1,1)=-athet(2,i,1,1)
+ bthet(1,i,-1,1)=bthet(1,i,1,1)
+ bthet(2,i,-1,1)=bthet(2,i,1,1)
+ enddo
+ do i=-ntyp,-1
+ a0thet(i)=a0thet(-i)
+ athet(1,i,-1,-1)=athet(1,-i,1,1)
+ athet(2,i,-1,-1)=-athet(2,-i,1,1)
+ bthet(1,i,-1,-1)=bthet(1,-i,1,1)
+ bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
+ athet(1,i,-1,1)=athet(1,-i,1,1)
+ athet(2,i,-1,1)=-athet(2,-i,1,1)
+ bthet(1,i,-1,1)=-bthet(1,-i,1,1)
+ bthet(2,i,-1,1)=bthet(2,-i,1,1)
+ athet(1,i,1,-1)=-athet(1,-i,1,1)
+ athet(2,i,1,-1)=athet(2,-i,1,1)
+ bthet(1,i,1,-1)=bthet(1,-i,1,1)
+ bthet(2,i,1,-1)=-bthet(2,-i,1,1)
+ theta0(i)=theta0(-i)
+ sig0(i)=sig0(-i)
+ sigc0(i)=sigc0(-i)
+ do j=0,3
+ polthet(j,i)=polthet(j,-i)
+ enddo
+ do j=1,3
+ gthet(j,i)=gthet(j,-i)
+ enddo
+ enddo
+c write (2,*) "End reading THETA_PDB"
+ close (ithep_pdb)
#endif
close(ithep)
#ifdef CRYST_SC
bsc(1,i)=0.0D0
read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
& ((blower(k,l,1),l=1,k),k=1,3)
+ censc(1,1,-i)=censc(1,1,i)
+ censc(2,1,-i)=censc(2,1,i)
+ censc(3,1,-i)=-censc(3,1,i)
do j=2,nlob(i)
read (irotam,*,end=112,err=112) bsc(j,i)
read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
& ((blower(k,l,j),l=1,k),k=1,3)
+ censc(1,j,-i)=censc(1,j,i)
+ censc(2,j,-i)=censc(2,j,i)
+ censc(3,j,-i)=-censc(3,j,i)
+C BSC is amplitude of Gaussian
enddo
do j=1,nlob(i)
do k=1,3
enddo
gaussc(k,l,j,i)=akl
gaussc(l,k,j,i)=akl
+ if (((k.eq.3).and.(l.ne.3))
+ & .or.((l.eq.3).and.(k.ne.3))) then
+ gaussc(k,l,j,-i)=-akl
+ gaussc(l,k,j,-i)=-akl
+ else
+ gaussc(k,l,j,-i)=akl
+ gaussc(l,k,j,-i)=akl
+ endif
enddo
enddo
enddo
enddo
endif
enddo
+C
+C Read the parameters of the probability distribution/energy expression
+C of the side chains.
+C
+ write (2,*) "Start reading ROTAM_PDB"
+ do i=1,ntyp
+ read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
+ if (i.eq.10) then
+ dsc_inv(i)=0.0D0
+ else
+ dsc_inv(i)=1.0D0/dsc(i)
+ endif
+ if (i.ne.10) then
+ do j=1,nlob(i)
+ do k=1,3
+ do l=1,3
+ blower(l,k,j)=0.0D0
+ enddo
+ enddo
+ enddo
+ bsc(1,i)=0.0D0
+ read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
+ & ((blower(k,l,1),l=1,k),k=1,3)
+ do j=2,nlob(i)
+ read (irotam_pdb,*,end=112,err=112) bsc(j,i)
+ read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
+ & ((blower(k,l,j),l=1,k),k=1,3)
+ enddo
+ do j=1,nlob(i)
+ do k=1,3
+ do l=1,k
+ akl=0.0D0
+ do m=1,3
+ akl=akl+blower(k,m,j)*blower(l,m,j)
+ enddo
+ gaussc(k,l,j,i)=akl
+ gaussc(l,k,j,i)=akl
+ enddo
+ enddo
+ enddo
+ endif
+ enddo
+ close (irotam_pdb)
+ write (2,*) "End reading ROTAM_PDB"
#endif
close(irotam)
+C
+C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
+C interaction energy of the Gly, Ala, and Pro prototypes.
+C
+ read (ifourier,*) nloctyp
+ SPLIT_FOURIERTOR = nloctyp.lt.0
+ nloctyp = iabs(nloctyp)
+#ifdef NEWCORR
+ read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
+ read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
+ itype2loc(ntyp1)=nloctyp
+ iloctyp(nloctyp)=ntyp1
+ do i=1,ntyp1
+ itype2loc(-i)=-itype2loc(i)
+ enddo
+#else
+ iloctyp(0)=10
+ iloctyp(1)=9
+ iloctyp(2)=20
+ iloctyp(3)=ntyp1
+#endif
+ do i=1,nloctyp
+ iloctyp(-i)=-iloctyp(i)
+ enddo
+c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
+c write (iout,*) "nloctyp",nloctyp,
+c & " iloctyp",(iloctyp(i),i=0,nloctyp)
+#ifdef NEWCORR
+ do i=0,nloctyp-1
+c write (iout,*) "NEWCORR",i
+ read (ifourier,*,end=115,err=115)
+ do ii=1,3
+ do j=1,2
+ read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
+ enddo
+ enddo
+c write (iout,*) "NEWCORR BNEW1"
+c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
+ do ii=1,3
+ do j=1,2
+ read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
+ enddo
+ enddo
+c write (iout,*) "NEWCORR BNEW2"
+c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
+ do kk=1,3
+ read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
+ read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
+ enddo
+c write (iout,*) "NEWCORR CCNEW"
+c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
+ do kk=1,3
+ read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
+ read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
+ enddo
+c write (iout,*) "NEWCORR DDNEW"
+c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
+ do ii=1,2
+ do jj=1,2
+ do kk=1,2
+ read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
+ enddo
+ enddo
+ enddo
+c write (iout,*) "NEWCORR EENEW1"
+c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
+ do ii=1,3
+ read (ifourier,*,end=115,err=115) e0new(ii,i)
+ enddo
+c write (iout,*) (e0new(ii,i),ii=1,3)
+ enddo
+c write (iout,*) "NEWCORR EENEW"
+ do i=0,nloctyp-1
+ do ii=1,3
+ ccnew(ii,1,i)=ccnew(ii,1,i)/2
+ ccnew(ii,2,i)=ccnew(ii,2,i)/2
+ ddnew(ii,1,i)=ddnew(ii,1,i)/2
+ ddnew(ii,2,i)=ddnew(ii,2,i)/2
+ enddo
+ enddo
+ do i=1,nloctyp-1
+ do ii=1,3
+ bnew1(ii,1,-i)= bnew1(ii,1,i)
+ bnew1(ii,2,-i)=-bnew1(ii,2,i)
+ bnew2(ii,1,-i)= bnew2(ii,1,i)
+ bnew2(ii,2,-i)=-bnew2(ii,2,i)
+ enddo
+ do ii=1,3
+c ccnew(ii,1,i)=ccnew(ii,1,i)/2
+c ccnew(ii,2,i)=ccnew(ii,2,i)/2
+c ddnew(ii,1,i)=ddnew(ii,1,i)/2
+c ddnew(ii,2,i)=ddnew(ii,2,i)/2
+ ccnew(ii,1,-i)=ccnew(ii,1,i)
+ ccnew(ii,2,-i)=-ccnew(ii,2,i)
+ ddnew(ii,1,-i)=ddnew(ii,1,i)
+ ddnew(ii,2,-i)=-ddnew(ii,2,i)
+ enddo
+ e0new(1,-i)= e0new(1,i)
+ e0new(2,-i)=-e0new(2,i)
+ e0new(3,-i)=-e0new(3,i)
+ do kk=1,2
+ eenew(kk,1,1,-i)= eenew(kk,1,1,i)
+ eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
+ eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
+ eenew(kk,2,2,-i)= eenew(kk,2,2,i)
+ enddo
+ enddo
+ if (lprint) then
+ write (iout,'(a)') "Coefficients of the multibody terms"
+ do i=-nloctyp+1,nloctyp-1
+ write (iout,*) "Type: ",onelet(iloctyp(i))
+ write (iout,*) "Coefficients of the expansion of B1"
+ do j=1,2
+ write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
+ enddo
+ write (iout,*) "Coefficients of the expansion of B2"
+ do j=1,2
+ write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
+ enddo
+ write (iout,*) "Coefficients of the expansion of C"
+ write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
+ write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
+ write (iout,*) "Coefficients of the expansion of D"
+ write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
+ write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
+ write (iout,*) "Coefficients of the expansion of E"
+ write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
+ do j=1,2
+ do k=1,2
+ write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
+ enddo
+ enddo
+ enddo
+ endif
+ IF (SPLIT_FOURIERTOR) THEN
+ do i=0,nloctyp-1
+c write (iout,*) "NEWCORR TOR",i
+ read (ifourier,*,end=115,err=115)
+ do ii=1,3
+ do j=1,2
+ read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i)
+ enddo
+ enddo
+c write (iout,*) "NEWCORR BNEW1 TOR"
+c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
+ do ii=1,3
+ do j=1,2
+ read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i)
+ enddo
+ enddo
+c write (iout,*) "NEWCORR BNEW2 TOR"
+c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
+ do kk=1,3
+ read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i)
+ read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i)
+ enddo
+c write (iout,*) "NEWCORR CCNEW TOR"
+c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
+ do kk=1,3
+ read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i)
+ read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i)
+ enddo
+c write (iout,*) "NEWCORR DDNEW TOR"
+c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
+ do ii=1,2
+ do jj=1,2
+ do kk=1,2
+ read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i)
+ enddo
+ enddo
+ enddo
+c write (iout,*) "NEWCORR EENEW1 TOR"
+c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
+ do ii=1,3
+ read (ifourier,*,end=115,err=115) e0newtor(ii,i)
+ enddo
+c write (iout,*) (e0newtor(ii,i),ii=1,3)
+ enddo
+c write (iout,*) "NEWCORR EENEW TOR"
+ do i=0,nloctyp-1
+ do ii=1,3
+ ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
+ ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
+ ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
+ ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
+ enddo
+ enddo
+ do i=1,nloctyp-1
+ do ii=1,3
+ bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
+ bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
+ bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
+ bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
+ enddo
+ do ii=1,3
+ ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
+ ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
+ ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
+ ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
+ enddo
+ e0newtor(1,-i)= e0newtor(1,i)
+ e0newtor(2,-i)=-e0newtor(2,i)
+ e0newtor(3,-i)=-e0newtor(3,i)
+ do kk=1,2
+ eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
+ eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
+ eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
+ eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
+ enddo
+ enddo
+ if (lprint) then
+ write (iout,'(a)')
+ & "Single-body coefficients of the torsional potentials"
+ do i=-nloctyp+1,nloctyp-1
+ write (iout,*) "Type: ",onelet(iloctyp(i))
+ write (iout,*) "Coefficients of the expansion of B1tor"
+ do j=1,2
+ write (iout,'(3hB1(,i1,1h),3f10.5)')
+ & j,(bnew1tor(k,j,i),k=1,3)
+ enddo
+ write (iout,*) "Coefficients of the expansion of B2tor"
+ do j=1,2
+ write (iout,'(3hB2(,i1,1h),3f10.5)')
+ & j,(bnew2tor(k,j,i),k=1,3)
+ enddo
+ write (iout,*) "Coefficients of the expansion of Ctor"
+ write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
+ write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
+ write (iout,*) "Coefficients of the expansion of Dtor"
+ write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
+ write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
+ write (iout,*) "Coefficients of the expansion of Etor"
+ write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
+ do j=1,2
+ do k=1,2
+ write (iout,'(1hE,2i1,2f10.5)')
+ & j,k,(eenewtor(l,j,k,i),l=1,2)
+ enddo
+ enddo
+ enddo
+ endif
+ ELSE
+ do i=-nloctyp+1,nloctyp-1
+ do ii=1,3
+ do j=1,2
+ bnew1tor(ii,j,i)=bnew1(ii,j,i)
+ enddo
+ enddo
+ do ii=1,3
+ do j=1,2
+ bnew2tor(ii,j,i)=bnew2(ii,j,i)
+ enddo
+ enddo
+ do ii=1,3
+ ccnewtor(ii,1,i)=ccnew(ii,1,i)
+ ccnewtor(ii,2,i)=ccnew(ii,2,i)
+ ddnewtor(ii,1,i)=ddnew(ii,1,i)
+ ddnewtor(ii,2,i)=ddnew(ii,2,i)
+ enddo
+ enddo
+ ENDIF !SPLIT_FOURIER_TOR
+#else
+ if (lprint)
+ & write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
+ do i=0,nloctyp-1
+ read (ifourier,*,end=115,err=115)
+ read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
+ if (lprint) then
+ write (iout,*) 'Type ',onelet(iloctyp(i))
+ write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
+ endif
+ if (i.gt.0) then
+ b(2,-i)= b(2,i)
+ b(3,-i)= b(3,i)
+ b(4,-i)=-b(4,i)
+ b(5,-i)=-b(5,i)
+ endif
+c B1(1,i) = b(3)
+c B1(2,i) = b(5)
+c B1(1,-i) = b(3)
+c B1(2,-i) = -b(5)
+c b1(1,i)=0.0d0
+c b1(2,i)=0.0d0
+c B1tilde(1,i) = b(3)
+c B1tilde(2,i) =-b(5)
+c B1tilde(1,-i) =-b(3)
+c B1tilde(2,-i) =b(5)
+c b1tilde(1,i)=0.0d0
+c b1tilde(2,i)=0.0d0
+c B2(1,i) = b(2)
+c B2(2,i) = b(4)
+c B2(1,-i) =b(2)
+c B2(2,-i) =-b(4)
+cc B1tilde(1,i) = b(3,i)
+cc B1tilde(2,i) =-b(5,i)
+C B1tilde(1,-i) =-b(3,i)
+C B1tilde(2,-i) =b(5,i)
+cc b1tilde(1,i)=0.0d0
+cc b1tilde(2,i)=0.0d0
+cc B2(1,i) = b(2,i)
+cc B2(2,i) = b(4,i)
+C B2(1,-i) =b(2,i)
+C B2(2,-i) =-b(4,i)
+
+c b2(1,i)=0.0d0
+c b2(2,i)=0.0d0
+ CCold(1,1,i)= b(7,i)
+ CCold(2,2,i)=-b(7,i)
+ CCold(2,1,i)= b(9,i)
+ CCold(1,2,i)= b(9,i)
+ CCold(1,1,-i)= b(7,i)
+ CCold(2,2,-i)=-b(7,i)
+ CCold(2,1,-i)=-b(9,i)
+ CCold(1,2,-i)=-b(9,i)
+c CC(1,1,i)=0.0d0
+c CC(2,2,i)=0.0d0
+c CC(2,1,i)=0.0d0
+c CC(1,2,i)=0.0d0
+c Ctilde(1,1,i)= CCold(1,1,i)
+c Ctilde(1,2,i)= CCold(1,2,i)
+c Ctilde(2,1,i)=-CCold(2,1,i)
+c Ctilde(2,2,i)=-CCold(2,2,i)
+
+c Ctilde(1,1,i)=0.0d0
+c Ctilde(1,2,i)=0.0d0
+c Ctilde(2,1,i)=0.0d0
+c Ctilde(2,2,i)=0.0d0
+ DDold(1,1,i)= b(6,i)
+ DDold(2,2,i)=-b(6,i)
+ DDold(2,1,i)= b(8,i)
+ DDold(1,2,i)= b(8,i)
+ DDold(1,1,-i)= b(6,i)
+ DDold(2,2,-i)=-b(6,i)
+ DDold(2,1,-i)=-b(8,i)
+ DDold(1,2,-i)=-b(8,i)
+c DD(1,1,i)=0.0d0
+c DD(2,2,i)=0.0d0
+c DD(2,1,i)=0.0d0
+c DD(1,2,i)=0.0d0
+c Dtilde(1,1,i)= DD(1,1,i)
+c Dtilde(1,2,i)= DD(1,2,i)
+c Dtilde(2,1,i)=-DD(2,1,i)
+c Dtilde(2,2,i)=-DD(2,2,i)
+
+c Dtilde(1,1,i)=0.0d0
+c Dtilde(1,2,i)=0.0d0
+c Dtilde(2,1,i)=0.0d0
+c Dtilde(2,2,i)=0.0d0
+ EEold(1,1,i)= b(10,i)+b(11,i)
+ EEold(2,2,i)=-b(10,i)+b(11,i)
+ EEold(2,1,i)= b(12,i)-b(13,i)
+ EEold(1,2,i)= b(12,i)+b(13,i)
+ EEold(1,1,-i)= b(10,i)+b(11,i)
+ EEold(2,2,-i)=-b(10,i)+b(11,i)
+ EEold(2,1,-i)=-b(12,i)+b(13,i)
+ EEold(1,2,-i)=-b(12,i)-b(13,i)
+ write(iout,*) "TU DOCHODZE"
+ print *,"JESTEM"
+c ee(1,1,i)=1.0d0
+c ee(2,2,i)=1.0d0
+c ee(2,1,i)=0.0d0
+c ee(1,2,i)=0.0d0
+c ee(2,1,i)=ee(1,2,i)
+ enddo
+ if (lprint) then
+ write (iout,*)
+ write (iout,*)
+ &"Coefficients of the cumulants (independent of valence angles)"
+ do i=-nloctyp+1,nloctyp-1
+ write (iout,*) 'Type ',onelet(iloctyp(i))
+ write (iout,*) 'B1'
+ write(iout,'(2f10.5)') B(3,i),B(5,i)
+ write (iout,*) 'B2'
+ write(iout,'(2f10.5)') B(2,i),B(4,i)
+ write (iout,*) 'CC'
+ do j=1,2
+ write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
+ enddo
+ write(iout,*) 'DD'
+ do j=1,2
+ write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
+ enddo
+ write(iout,*) 'EE'
+ do j=1,2
+ write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
+ enddo
+ enddo
+ endif
+#endif
#ifdef CRYST_TOR
C
C
C Read torsional parameters
C
+ IF (TOR_MODE.eq.0) THEN
+
read (itorp,*,end=113,err=113) ntortyp
read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
-c write (iout,*) 'ntortyp',ntortyp
- do i=1,ntortyp
- do j=1,ntortyp
- read (itorp,*,end=113,err=113) nterm(i,j),nlor(i,j)
+ do i=1,ntyp1
+ itype2loc(i)=itortyp(i)
+ enddo
+ do i=1,ntyp1
+ itype2loc(-i)=-itype2loc(i)
+ enddo
+ itortyp(ntyp1)=ntortyp
+ do iblock=1,2
+ do i=-ntyp,-1
+ itortyp(i)=-itortyp(-i)
+ enddo
+ write (iout,*) 'ntortyp',ntortyp
+ do i=0,ntortyp-1
+ do j=-ntortyp+1,ntortyp-1
+ read (itorp,*,end=113,err=113) nterm(i,j,iblock),
+ & nlor(i,j,iblock)
+ nterm(-i,-j,iblock)=nterm(i,j,iblock)
+ nlor(-i,-j,iblock)=nlor(i,j,iblock)
v0ij=0.0d0
si=-1.0d0
- do k=1,nterm(i,j)
- read (itorp,*,end=113,err=113) kk,v1(k,i,j),v2(k,i,j)
- v0ij=v0ij+si*v1(k,i,j)
+ do k=1,nterm(i,j,iblock)
+ read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
+ & v2(k,i,j,iblock)
+ v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
+ v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
+ v0ij=v0ij+si*v1(k,i,j,iblock)
si=-si
+c write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
+c write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
+c &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
enddo
- do k=1,nlor(i,j)
+ do k=1,nlor(i,j,iblock)
read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
- & vlor2(k,i,j),vlor3(k,i,j)
+ & vlor2(k,i,j),vlor3(k,i,j)
v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
enddo
- v0(i,j)=v0ij
+ v0(i,j,iblock)=v0ij
+ v0(-i,-j,iblock)=v0ij
enddo
enddo
+ enddo
close (itorp)
if (lprint) then
- write (iout,'(/a/)') 'Torsional constants:'
- do i=1,ntortyp
- do j=1,ntortyp
+ write (iout,'(/a/)') 'Parameters of the SCCOR potentials:'
+ do iblock=1,2
+ do i=0,ntortyp-1
+ do j=-ntortyp+1,ntortyp-1
write (iout,*) 'ityp',i,' jtyp',j
write (iout,*) 'Fourier constants'
- do k=1,nterm(i,j)
- write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
+ do k=1,nterm(i,j,iblock)
+ write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
+ & v2(k,i,j,iblock)
enddo
write (iout,*) 'Lorenz constants'
- do k=1,nlor(i,j)
- write (iout,'(3(1pe15.5))')
+ do k=1,nlor(i,j,iblock)
+ write (iout,'(3(1pe15.5))')
& vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
enddo
enddo
enddo
+ enddo
endif
+
C
C 6/23/01 Read parameters for double torsionals
C
- do i=1,ntortyp
- do j=1,ntortyp
- do k=1,ntortyp
+ do iblock=1,2
+ do i=0,ntortyp-1
+ do j=-ntortyp+1,ntortyp-1
+ do k=-ntortyp+1,ntortyp-1
read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
- if (t1.ne.onelett(i) .or. t2.ne.onelett(j)
- & .or. t3.ne.onelett(k)) then
+c write (iout,*) "OK onelett",
+c & i,j,k,t1,t2,t3
+
+ if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
+ & .or. t3.ne.toronelet(k)) then
write (iout,*) "Error in double torsional parameter file",
& i,j,k,t1,t2,t3
#ifdef MPI
#endif
stop "Error in double torsional parameter file"
endif
- read (itordp,*,end=114,err=114) ntermd_1(i,j,k),
- & ntermd_2(i,j,k)
- read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k),l=1,
- & ntermd_1(i,j,k))
- read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k),l=1,
- & ntermd_1(i,j,k))
- read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k),l=1,
- & ntermd_1(i,j,k))
- read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k),l=1,
- & ntermd_1(i,j,k))
- read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k),
- & v2c(m,l,i,j,k),v2s(l,m,i,j,k),v2s(m,l,i,j,k),
- & m=1,l-1),l=1,ntermd_2(i,j,k))
- enddo
- enddo
- enddo
+ read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
+ & ntermd_2(i,j,k,iblock)
+ ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
+ ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
+ read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
+ & ntermd_1(i,j,k,iblock))
+ read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
+ & ntermd_1(i,j,k,iblock))
+ read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
+ & ntermd_1(i,j,k,iblock))
+ read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
+ & ntermd_1(i,j,k,iblock))
+C Martix of D parameters for one dimesional foureir series
+ do l=1,ntermd_1(i,j,k,iblock)
+ v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
+ v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
+ v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
+ v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
+c write(iout,*) "whcodze" ,
+c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
+ enddo
+ read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
+ & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
+ & v2s(m,l,i,j,k,iblock),
+ & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
+C Martix of D parameters for two dimesional fourier series
+ do l=1,ntermd_2(i,j,k,iblock)
+ do m=1,l-1
+ v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
+ v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
+ v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
+ v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
+ enddo!m
+ enddo!l
+ enddo!k
+ enddo!j
+ enddo!i
+ enddo!iblock
if (lprint) then
- write (iout,*)
+ write (iout,*)
write (iout,*) 'Constants for double torsionals'
- do i=1,ntortyp
- do j=1,ntortyp
- do k=1,ntortyp
+ do iblock=1,2
+ do i=0,ntortyp-1
+ do j=-ntortyp+1,ntortyp-1
+ do k=-ntortyp+1,ntortyp-1
write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
- & ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k)
+ & ' nsingle',ntermd_1(i,j,k,iblock),
+ & ' ndouble',ntermd_2(i,j,k,iblock)
write (iout,*)
write (iout,*) 'Single angles:'
- do l=1,ntermd_1(i,j,k)
- write (iout,'(i5,2f10.5,5x,2f10.5)') l,
- & v1c(1,l,i,j,k),v1s(1,l,i,j,k),
- & v1c(2,l,i,j,k),v1s(2,l,i,j,k)
+ do l=1,ntermd_1(i,j,k,iblock)
+ write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
+ & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
+ & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
+ & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
enddo
write (iout,*)
write (iout,*) 'Pairs of angles:'
- write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
- do l=1,ntermd_2(i,j,k)
- write (iout,'(i5,20f10.5)')
- & l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k))
+ write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
+ do l=1,ntermd_2(i,j,k,iblock)
+ write (iout,'(i5,20f10.5)')
+ & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
enddo
write (iout,*)
- write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
- do l=1,ntermd_2(i,j,k)
- write (iout,'(i5,20f10.5)')
- & l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k))
+ write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
+ do l=1,ntermd_2(i,j,k,iblock)
+ write (iout,'(i5,20f10.5)')
+ & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
+ & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
enddo
write (iout,*)
enddo
enddo
enddo
+ enddo
endif
+
+ ELSE IF (TOR_MODE.eq.1) THEN
+
+C read valence-torsional parameters
+ read (itorp,*,end=121,err=121) ntortyp
+ nkcctyp=ntortyp
+ write (iout,*) "Valence-torsional parameters read in ntortyp",
+ & ntortyp
+ read (itorp,*,end=121,err=121) (itortyp(i),i=1,ntyp)
+ write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
+#ifndef NEWCORR
+ do i=1,ntyp1
+ itype2loc(i)=itortyp(i)
+ enddo
#endif
-C
-C 5/21/07 (AL) Read coefficients of the backbone-local sidechain-local
-C correlation energies.
-C
- read (isccor,*,end=119,err=119) nterm_sccor
- do i=1,20
- do j=1,20
- read (isccor,'(a)')
- do k=1,nterm_sccor
- read (isccor,*,end=119,err=119) kk,v1sccor(k,i,j),
- & v2sccor(k,i,j)
+ do i=-ntyp,-1
+ itortyp(i)=-itortyp(-i)
+ enddo
+ do i=-ntortyp+1,ntortyp-1
+ do j=-ntortyp+1,ntortyp-1
+C first we read the cos and sin gamma parameters
+ read (itorp,'(13x,a)',end=121,err=121) string
+ write (iout,*) i,j,string
+ read (itorp,*,end=121,err=121)
+ & nterm_kcc(j,i),nterm_kcc_Tb(j,i)
+C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
+ do k=1,nterm_kcc(j,i)
+ do l=1,nterm_kcc_Tb(j,i)
+ do ll=1,nterm_kcc_Tb(j,i)
+ read (itorp,*,end=121,err=121) ii,jj,kk,
+ & v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
+ enddo
+ enddo
enddo
enddo
enddo
- close (isccor)
- if (lprint) then
- write (iout,'(/a/)') 'Torsional constants of SCCORR:'
- do i=1,20
- do j=1,20
- write (iout,*) 'ityp',i,' jtyp',j
- do k=1,nterm_sccor
- write (iout,'(2(1pe15.5))') v1sccor(k,i,j),v2sccor(k,i,j)
+ ELSE
+#ifdef NEWCORR
+c AL 4/8/16: Calculate coefficients from one-body parameters
+ ntortyp=nloctyp
+ do i=-ntyp1,ntyp1
+ itortyp(i)=itype2loc(i)
+ enddo
+ write (iout,*)
+ &"Val-tor parameters calculated from cumulant coefficients ntortyp"
+ & ,ntortyp
+ do i=-ntortyp+1,ntortyp-1
+ do j=-ntortyp+1,ntortyp-1
+ nterm_kcc(j,i)=2
+ nterm_kcc_Tb(j,i)=3
+ do k=1,nterm_kcc_Tb(j,i)
+ do l=1,nterm_kcc_Tb(j,i)
+ v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)
+ & +bnew1tor(k,2,i)*bnew2tor(l,2,j)
+ v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)
+ & +bnew1tor(k,2,i)*bnew2tor(l,1,j)
enddo
enddo
+ do k=1,nterm_kcc_Tb(j,i)
+ do l=1,nterm_kcc_Tb(j,i)
+#ifdef CORRCD
+ v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j)
+ & -ccnewtor(k,2,i)*ddnewtor(l,2,j))
+ v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j)
+ & +ccnewtor(k,1,i)*ddnewtor(l,2,j))
+#else
+ v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j)
+ & -ccnewtor(k,2,i)*ddnewtor(l,2,j))
+ v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j)
+ & +ccnewtor(k,1,i)*ddnewtor(l,2,j))
+#endif
+ enddo
+ enddo
+c f(theta,gamma)=-(b21(theta)*b11(theta)+b12(theta)*b22(theta))*cos(gamma)-(b22(theta)*b11(theta)+b21(theta)*b12(theta))*sin(gamma)+(c11(theta)*d11(theta)-c12(theta)*d12(theta))*cos(2*gamma)+(c12(theta)*d11(theta)+c11(theta)*d12(theta))*sin(2*gamma)
+ enddo
+ enddo
+#else
+ write (iout,*) "TOR_MODE>1 only with NEWCORR"
+ stop
+#endif
+ ENDIF ! TOR_MODE
+
+ if (tor_mode.gt.0 .and. lprint) then
+c Print valence-torsional parameters
+ write (iout,'(a)')
+ & "Parameters of the valence-torsional potentials"
+ do i=-ntortyp+1,ntortyp-1
+ do j=-ntortyp+1,ntortyp-1
+ write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
+ write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
+ do k=1,nterm_kcc(j,i)
+ do l=1,nterm_kcc_Tb(j,i)
+ do ll=1,nterm_kcc_Tb(j,i)
+ write (iout,'(3i5,2f15.4)')
+ & k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
+ enddo
+ enddo
+ enddo
+ enddo
enddo
endif
+
+#endif
+C Read of Side-chain backbone correlation parameters
+C Modified 11 May 2012 by Adasko
+CCC
C
-C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
-C interaction energy of the Gly, Ala, and Pro prototypes.
-C
- if (lprint) then
- write (iout,*)
- write (iout,*) "Coefficients of the cumulants"
- endif
- read (ifourier,*) nloctyp
- do i=1,nloctyp
- read (ifourier,*,end=115,err=115)
- read (ifourier,*,end=115,err=115) (b(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)
- endif
- B1(1,i) = b(3)
- B1(2,i) = b(5)
-c b1(1,i)=0.0d0
-c b1(2,i)=0.0d0
- B1tilde(1,i) = b(3)
- B1tilde(2,i) =-b(5)
-c b1tilde(1,i)=0.0d0
-c b1tilde(2,i)=0.0d0
- B2(1,i) = b(2)
- B2(2,i) = b(4)
-c b2(1,i)=0.0d0
-c 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)
-c CC(1,1,i)=0.0d0
-c CC(2,2,i)=0.0d0
-c CC(2,1,i)=0.0d0
-c 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)
-c Ctilde(1,1,i)=0.0d0
-c Ctilde(1,2,i)=0.0d0
-c Ctilde(2,1,i)=0.0d0
-c 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)
-c DD(1,1,i)=0.0d0
-c DD(2,2,i)=0.0d0
-c DD(2,1,i)=0.0d0
-c 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)
-c Dtilde(1,1,i)=0.0d0
-c Dtilde(1,2,i)=0.0d0
-c Dtilde(2,1,i)=0.0d0
-c 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)
-c ee(1,1,i)=1.0d0
-c ee(2,2,i)=1.0d0
-c ee(2,1,i)=0.0d0
-c ee(1,2,i)=0.0d0
-c ee(2,1,i)=ee(1,2,i)
+ read (isccor,*,end=119,err=119) nsccortyp
+#ifdef SCCORPDB
+ read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
+ do i=-ntyp,-1
+ isccortyp(i)=-isccortyp(-i)
enddo
- if (lprint) then
- do i=1,nloctyp
- write (iout,*) 'Type',i
- write (iout,*) 'B1'
- write(iout,*) B1(1,i),B1(2,i)
- write (iout,*) 'B2'
- write(iout,*) B2(1,i),B2(2,i)
- write (iout,*) 'CC'
- do j=1,2
- write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
- enddo
- write(iout,*) 'DD'
- do j=1,2
- write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
+ iscprol=isccortyp(20)
+c write (iout,*) 'ntortyp',ntortyp
+ maxinter=3
+cc maxinter is maximum interaction sites
+ do l=1,maxinter
+ do i=1,nsccortyp
+ do j=1,nsccortyp
+ read (isccor,*,end=119,err=119)
+ &nterm_sccor(i,j),nlor_sccor(i,j)
+ v0ijsccor=0.0d0
+ v0ijsccor1=0.0d0
+ v0ijsccor2=0.0d0
+ v0ijsccor3=0.0d0
+ si=-1.0d0
+ nterm_sccor(-i,j)=nterm_sccor(i,j)
+ nterm_sccor(-i,-j)=nterm_sccor(i,j)
+ nterm_sccor(i,-j)=nterm_sccor(i,j)
+ do k=1,nterm_sccor(i,j)
+ read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
+ & ,v2sccor(k,l,i,j)
+ if (j.eq.iscprol) then
+ if (i.eq.isccortyp(10)) then
+ v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
+ v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
+ else
+ v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
+ & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
+ v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
+ & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
+ v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
+ v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
+ v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
+ v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
+ endif
+ else
+ if (i.eq.isccortyp(10)) then
+ v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
+ v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
+ else
+ if (j.eq.isccortyp(10)) then
+ v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
+ v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
+ else
+ v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
+ v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
+ v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
+ v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
+ v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
+ v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
+ endif
+ endif
+ endif
+ v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
+ v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
+ v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
+ v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
+ si=-si
+ enddo
+ do k=1,nlor_sccor(i,j)
+ read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
+ & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
+ v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
+ &(1+vlor3sccor(k,i,j)**2)
+ enddo
+ v0sccor(l,i,j)=v0ijsccor
+ v0sccor(l,-i,j)=v0ijsccor1
+ v0sccor(l,i,-j)=v0ijsccor2
+ v0sccor(l,-i,-j)=v0ijsccor3
enddo
- write(iout,*) 'EE'
- do j=1,2
- write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
+ enddo
+ enddo
+ close (isccor)
+#else
+ read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
+c write (iout,*) 'ntortyp',ntortyp
+ maxinter=3
+cc maxinter is maximum interaction sites
+ do l=1,maxinter
+ do i=1,nsccortyp
+ do j=1,nsccortyp
+ read (isccor,*,end=119,err=119)
+ & nterm_sccor(i,j),nlor_sccor(i,j)
+ v0ijsccor=0.0d0
+ si=-1.0d0
+
+ do k=1,nterm_sccor(i,j)
+ read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
+ & ,v2sccor(k,l,i,j)
+ v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
+ si=-si
+ enddo
+ do k=1,nlor_sccor(i,j)
+ read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
+ & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
+ v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
+ &(1+vlor3sccor(k,i,j)**2)
+ enddo
+ v0sccor(l,i,j)=v0ijsccor
enddo
enddo
+ enddo
+ close (isccor)
+
+#endif
+ if (lprint) then
+ write (iout,'(/a/)') 'Torsional constants:'
+ do l=1,maxinter
+ do i=1,nsccortyp
+ do j=1,nsccortyp
+ write (iout,*) 'ityp',i,' jtyp',j
+ write (iout,*) 'Fourier constants'
+ do k=1,nterm_sccor(i,j)
+ write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
+ enddo
+ write (iout,*) 'Lorenz constants'
+ do k=1,nlor_sccor(i,j)
+ write (iout,'(3(1pe15.5))')
+ & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
+ enddo
+ enddo
+ enddo
+ enddo
endif
+
C
C Read electrostatic-interaction parameters
C
bpp (i,j)=-2.0D0*epp(i,j)*rri
ael6(i,j)=elpp6(i,j)*4.2D0**6
ael3(i,j)=elpp3(i,j)*4.2D0**3
+c lprint=.true.
if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
& ael6(i,j),ael3(i,j)
+c lprint=.false.
enddo
enddo
C
& ', exponents are ',expon,2*expon
goto (10,20,30,30,40) ipot
C----------------------- LJ potential ---------------------------------
- 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
+ 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
& (sigma0(i),i=1,ntyp)
if (lprint) then
write (iout,'(/a/)') 'Parameters of the LJ potential:'
endif
goto 50
C----------------------- LJK potential --------------------------------
- 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
+ 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
& (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
if (lprint) then
write (iout,'(/a/)') 'Parameters of the LJK potential:'
endif
goto 50
C---------------------- GB or BP potential -----------------------------
- 30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
- & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
- & (alp(i),i=1,ntyp)
+ 30 do i=1,ntyp
+ read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
+ enddo
+ read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
+ 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)
+C now we start reading lipid
+ do i=1,ntyp
+ read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
+ if (lprint) write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp)
+
+C print *,"WARNING!!"
+C do j=1,ntyp
+C epslip(i,j)=epslip(i,j)+0.05d0
+C enddo
+ enddo
+ if (lprint) write(iout,*) epslip(1,1),"OK?"
C For the GB potential convert sigma'**2 into chi'
if (ipot.eq.4) then
- do i=1,ntyp
- chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
+ do i=1,ntyp
+ chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
enddo
endif
if (lprint) then
- write (iout,'(/a/)') 'Parameters of the BP potential:'
- write (iout,'(a/)') 'The epsilon array:'
- call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
- write (iout,'(/a)') 'One-body parameters:'
- write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
+ write (iout,'(/a/)') 'Parameters of the BP potential:'
+ write (iout,'(a/)') 'The epsilon array:'
+ call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
+ 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),sigma0(i),sigii(i),
& chip(i),alp(i),i=1,ntyp)
endif
goto 50
C--------------------- GBV potential -----------------------------------
- 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
- & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
- & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
+c 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
+c & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
+c & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
+ 40 do i=1,ntyp
+ read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
+ enddo
+ read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
+ read (isidep,*,end=117,err=117)(rr0(i),i=1,ntyp)
+ 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=1161,err=1161)(epslip(i,j),j=i,ntyp)
+ if (lprint) write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp)
+ enddo
+ do i=1,ntyp
+ chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
+ enddo
if (lprint) then
- write (iout,'(/a/)') 'Parameters of the GBV potential:'
- write (iout,'(a/)') 'The epsilon array:'
- call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
- write (iout,'(/a)') 'One-body parameters:'
- write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
+ write (iout,'(/a/)') 'Parameters of the GBV potential:'
+ write (iout,'(a/)') 'The epsilon array:'
+ call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
+ 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),sigma0(i),rr0(i),
& sigii(i),chip(i),alp(i),i=1,ntyp)
endif
+C now we start reading lipid
50 continue
close (isidep)
C-----------------------------------------------------------------------
C Calculate the "working" parameters of SC interactions.
do i=2,ntyp
do j=1,i-1
- eps(i,j)=eps(j,i)
+ 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)
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),restyp(j),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
+#ifdef TUBE
+ write(iout,*) "tube param"
+ read(itube,*) epspeptube,sigmapeptube
+ 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
+ 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
+#endif
+
#ifdef OLDSCP
C
C Define the SC-p interaction constants (hard-coded; old style)
C
- do i=1,20
+ do i=1,ntyp
C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
C helix formation)
c aad(i,1)=0.3D0*4.0D0**12
bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
enddo
-
+c lprint=.true.
if (lprint) then
write (iout,*) "Parameters of SC-p interactions:"
- do i=1,20
+ do i=1,ntyp
write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
& eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
enddo
endif
+c lprint=.false.
#endif
C
C Define the constants of the disulfide bridge
C
- ebr=-5.50D0
+C ebr=-12.00D0
c
c Old arbitrary potential - commented out.
c
c energy surface of diethyl disulfide.
c A. Liwo and U. Kozlowska, 11/24/03
c
- D0CM = 3.78d0
- AKCM = 15.1d0
- AKTH = 11.0d0
- AKCT = 12.0d0
- V1SS =-1.08d0
- V2SS = 7.61d0
- V3SS = 13.7d0
+C D0CM = 3.78d0
+C AKCM = 15.1d0
+C AKTH = 11.0d0
+C AKCT = 12.0d0
+C V1SS =-1.08d0
+C V2SS = 7.61d0
+C V3SS = 13.7d0
c akcm=0.0d0
c akth=0.0d0
c akct=0.0d0
c v2ss=0.0d0
c v3ss=0.0d0
- if(me.eq.king) then
- write (iout,'(/a)') "Disulfide bridge parameters:"
- write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
- write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
- write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
- write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
- & ' v3ss:',v3ss
+C if(me.eq.king) then
+C write (iout,'(/a)') "Disulfide bridge parameters:"
+C write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
+C write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
+C write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
+C write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
+C & ' v3ss:',v3ss
+C endif
+ if (shield_mode.gt.0) then
+C VSolvSphere the volume of solving sphere
+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
+c write (iout,*) "rpp(1,1)",rpp(1,1)," pi",pi
+ 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,"VSolvSphere_div",
+ & 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)
+ enddo
+ buff_shield=1.0d0
endif
return
111 write (iout,*) "Error reading bending energy parameters."
goto 999
116 write (iout,*) "Error reading electrostatic energy parameters."
goto 999
+ 1161 write (iout,*) "Error reading electrostatic energy parameters.Lip"
+ goto 999
117 write (iout,*) "Error reading side chain interaction parameters."
goto 999
118 write (iout,*) "Error reading SCp interaction parameters."
goto 999
119 write (iout,*) "Error reading SCCOR parameters"
+ go to 999
+ 121 write (iout,*) "Error in Czybyshev parameters"
999 continue
#ifdef MPI
call MPI_Finalize(Ierror)
#else
call getenv(var,val)
#endif
-
+C set the variables used for shielding effect
+C if (shield_mode.gt.0) then
+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
+C VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
+C VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
+C long axis of side chain
+C do i=1,ntyp
+C long_r_sidechain(i)=vbldsc0(1,i)
+C short_r_sidechain(i)=sigma0(i)
+C enddo
+C lets set the buffor value
+C buff_shield=1.0d0
+C endif
return
end
projects / unres.git / blobdiff