enddo
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)
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 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
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
enddo
call flush(iout)
endif
+ 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
+ write (2,*) "End reading THETA_PDB"
+ close (ithep_pdb)
#endif
close(ithep)
#ifdef CRYST_SC
C added by Urszula Kozlowska 07/11/2007
C
do i=1,ntyp
-CC TU JEST ZLE musibyc ntyp
read (irotam,*,end=112,err=112)
if (i.eq.10) then
read (irotam,*,end=112,err=112)
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)
si=-si
enddo
do k=1,nlor_sccor(i,j)
- read (isccor,*,end=113,err=113) kk,vlor1sccor(k,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
close (isccor)
#else
- read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
+ 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=113,err=113)
+ 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=113,err=113) kk,v1sccor(k,l,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=113,err=113) kk,vlor1sccor(k,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)
& ', 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:'
goto 50
C---------------------- GB or BP potential -----------------------------
30 do i=1,ntyp
- read (isidep,*,end=116,err=116)(eps(i,j),j=i,ntyp)
+ read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
enddo
- read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
- read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
- read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
- read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
+ 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 For the GB potential convert sigma'**2 into chi'
if (ipot.eq.4) then
do 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),
+ 40 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),(sigii(i),i=1,ntyp),
& (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
if (lprint) then
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
- endif
+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
return
111 write (iout,*) "Error reading bending energy parameters."
goto 999