& ntheterm3,nsingle,ndouble
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
+c VAR: aa0thet is variable describing the average value of Foureir
+c VAR: expansion series
+ read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
+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
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
do i=1,nthetyp+1
do j=1,nthetyp+1
do k=1,nthetyp+1
- write (iout,'(//4a)')
- & 'Type ',onelett(i),onelett(j),onelett(k)
+ write (iout,'(//4a)')
+ & 'Type ',onelett(i),onelett(j),onelett(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]))')
+ 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
write (iout,*)
write (iout,*) 'Constants for double torsionals'
do iblock=1,2
- do i=1,ntortyp
- do j=-ntortyp,ntortyp
- do k=-ntortyp,ntortyp
+ 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,iblock),
& ' ndouble',ntermd_2(i,j,k,iblock)
C Modified 11 May 2012 by Adasko
CCC
C
- read (isccor,*,end=113,err=113) nsccortyp
+ read (isccor,*,end=119,err=119) nsccortyp
#ifdef SCCORPDB
- read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
+ read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
do i=-ntyp,-1
isccortyp(i)=-isccortyp(-i)
enddo
do l=1,maxinter
do i=1,nsccortyp
do j=1,nsccortyp
- read (isccor,*,end=113,err=113) nterm_sccor(i,j),nlor_sccor(i,j)
+ read (isccor,*,end=119,err=119) nterm_sccor(i,j),nlor_sccor(i,j)
v0ijsccor=0.0d0
v0ijsccor1=0.0d0
v0ijsccor2=0.0d0
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=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)
if (j.eq.iscprol) then
if (i.eq.isccortyp(10)) then
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)
enddo
- v0sccor(i,j)=v0ijsccor
+ v0sccor(i,j,iblock)=v0ijsccor
enddo
enddo
enddo
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=116,err=116)(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)
+
+c 30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
+c & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
+c & (alp(i),i=1,ntyp)
C For the GB potential convert sigma'**2 into chi'
if (ipot.eq.4) then
do i=1,ntyp
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
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
projects / unres.git / blobdiff