From: Adam Kazimierz Sieradzan Date: Sat, 31 Aug 2013 12:22:46 +0000 (-0400) Subject: Clean-up in parmread for src-MD and src-MD-M X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=commitdiff_plain;h=3e2116f1c3c9f9db1f8f34ea327e064df9345d67 Clean-up in parmread for src-MD and src-MD-M pre-preparation for removal of src-MD and leaving out only the multichain version --- diff --git a/.gitignore b/.gitignore index 77eb558..b25ea14 100644 --- a/.gitignore +++ b/.gitignore @@ -10,6 +10,7 @@ # ignore build dir build/ +build2/ # ignored dirs form adasko gradcheck/ diff --git a/source/unres/src_MD-M/COMMON.LOCAL b/source/unres/src_MD-M/COMMON.LOCAL index 07ac66d..9f0627b 100644 --- a/source/unres/src_MD-M/COMMON.LOCAL +++ b/source/unres/src_MD-M/COMMON.LOCAL @@ -9,7 +9,7 @@ C Parameters of the virtual-bond-angle probability distribution C Parameters of the side-chain probability distribution common /sclocal/ dsc(ntyp1),dsc_inv(ntyp1),bsc(maxlob,ntyp), & censc(3,maxlob,-ntyp:ntyp),gaussc(3,3,maxlob,-ntyp:ntyp), - &d sc0(ntyp1), + & dsc0(ntyp1), & nlob(ntyp1) C Parameters of ab initio-derived potential of virtual-bond-angle bending integer nthetyp,ntheterm,ntheterm2,ntheterm3,nsingle,ndouble, diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F index 44998e3..154fcf4 100644 --- a/source/unres/src_MD-M/parmread.F +++ b/source/unres/src_MD-M/parmread.F @@ -243,12 +243,21 @@ C 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) @@ -269,7 +278,8 @@ C 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 @@ -285,6 +295,7 @@ 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 @@ -300,7 +311,7 @@ c enddo 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 @@ -377,6 +388,51 @@ C enddo call flush(iout) endif + write (2,*) "Start reading THETA_PDB" + do i=1,ntyp + 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 + 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 @@ -470,7 +526,6 @@ C Read scrot parameters for potentials determined from all-atom AM1 calculations 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) @@ -480,6 +535,48 @@ CC TU JEST ZLE musibyc ntyp enddo endif enddo +C +C Read the parameters of the probability distribution/energy expression +C of the side chains. +C + 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) #endif close(irotam) @@ -734,7 +831,7 @@ cc maxinter is maximum interaction sites 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) @@ -748,26 +845,26 @@ cc maxinter is maximum interaction sites 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) diff --git a/source/unres/src_MD/COMMON.LOCAL b/source/unres/src_MD/COMMON.LOCAL index f8b5ab3..9f0627b 100644 --- a/source/unres/src_MD/COMMON.LOCAL +++ b/source/unres/src_MD/COMMON.LOCAL @@ -3,27 +3,33 @@ integer nlob C Parameters of the virtual-bond-angle probability distribution common /thetas/ a0thet(-ntyp:ntyp),athet(2,-ntyp:ntyp,-1:1,-1:1), - & bthet(2,-ntyp:ntyp,-1:1,-1:1), polthet(0:3,-ntyp:ntyp), - & gthet(3,-ntyp:ntyp),theta0(-ntyp:ntyp),sig0(-ntyp:ntyp), + & bthet(2,-ntyp:ntyp,-1:1,-1:1),polthet(0:3,-ntyp:ntyp), + & gthet(3,-ntyp:ntyp),theta0(-ntyp:ntyp),sig0(-ntyp:ntyp), & sigc0(-ntyp:ntyp) C Parameters of the side-chain probability distribution common /sclocal/ dsc(ntyp1),dsc_inv(ntyp1),bsc(maxlob,ntyp), & censc(3,maxlob,-ntyp:ntyp),gaussc(3,3,maxlob,-ntyp:ntyp), - & dsc0(ntyp1), + & dsc0(ntyp1), & nlob(ntyp1) C Parameters of ab initio-derived potential of virtual-bond-angle bending integer nthetyp,ntheterm,ntheterm2,ntheterm3,nsingle,ndouble, - & ithetyp(ntyp1),nntheterm - double precision aa0thet(maxthetyp1,maxthetyp1,maxthetyp1), - & aathet(maxtheterm,maxthetyp1,maxthetyp1,maxthetyp1), - & bbthet(maxsingle,maxtheterm2,maxthetyp1,maxthetyp1,maxthetyp1), - & ccthet(maxsingle,maxtheterm2,maxthetyp1,maxthetyp1,maxthetyp1), - & ddthet(maxsingle,maxtheterm2,maxthetyp1,maxthetyp1,maxthetyp1), - & eethet(maxsingle,maxtheterm2,maxthetyp1,maxthetyp1,maxthetyp1), - & ffthet(maxdouble,maxdouble,maxtheterm3,maxthetyp1,maxthetyp1, - & maxthetyp1), - & ggthet(maxdouble,maxdouble,maxtheterm3,maxthetyp1,maxthetyp1, - & maxthetyp1) + & ithetyp(-ntyp1:ntyp1),nntheterm + double precision aa0thet(-maxthetyp1:maxthetyp1, + &-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2), + & aathet(maxtheterm,-maxthetyp1:maxthetyp1, + &-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2), + & bbthet(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1, + &-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2), + & ccthet(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1, + &-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2), + & ddthet(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1, + &-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2), + & eethet(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1, + &-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2), + & ffthet(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1, + &-maxthetyp1:maxthetyp1, -maxthetyp1:maxthetyp1,2), + & ggthet(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1, + &-maxthetyp1:maxthetyp1, -maxthetyp1:maxthetyp1,2) common /theta_abinitio/aa0thet,aathet,bbthet,ccthet,ddthet,eethet, & ffthet, & ggthet,ithetyp,nthetyp,ntheterm,ntheterm2,ntheterm3,nsingle, @@ -31,11 +37,10 @@ C Parameters of ab initio-derived potential of virtual-bond-angle bending C Virtual-bond lenghts double precision vbl,vblinv,vblinv2,vbl_cis,vbl0,vbld_inv integer loc_start,loc_end,ithet_start,ithet_end,iphi_start, - & iphi_end,iphid_start,iphid_end,itau_start,itau_end,ibond_start, - & ibond_end, + & iphi_end,iphid_start,iphid_end,ibond_start,ibond_end, & ibondp_start,ibondp_end,ivec_start,ivec_end,iset_start,iset_end, & iturn3_start,iturn3_end,iturn4_start,iturn4_end,iint_start, - & iint_end,iphi1_start,iphi1_end, + & iint_end,iphi1_start,iphi1_end,itau_start,itau_end, & ibond_displ(0:max_fg_procs-1),ibond_count(0:max_fg_procs-1), & ithet_displ(0:max_fg_procs-1),ithet_count(0:max_fg_procs-1), & iphi_displ(0:max_fg_procs-1),iphi_count(0:max_fg_procs-1), @@ -45,12 +50,11 @@ C Virtual-bond lenghts & iint_count(0:max_fg_procs-1),iint_displ(0:max_fg_procs-1) common /peptbond/ vbl,vblinv,vblinv2,vbl_cis,vbl0 common /indices/ loc_start,loc_end,ithet_start,ithet_end, - & iphi_start,iphi_end,iphid_start,iphid_end,itau_start,itau_end, - & ibond_start,ibond_end, + & iphi_start,iphi_end,iphid_start,iphid_end,ibond_start,ibond_end, & ibondp_start,ibondp_end,ivec_start,ivec_end,iset_start,iset_end, & iturn3_start,iturn3_end,iturn4_start,iturn4_end,iint_start, & iint_end,iphi1_start,iphi1_end,iint_count,iint_displ,ivec_displ, - & ivec_count,iset_displ, + & ivec_count,iset_displ,itau_start,itau_end, & iset_count,ibond_displ,ibond_count,ithet_displ,ithet_count, & iphi_displ,iphi_count,iphi1_displ,iphi1_count C Inverses of the actual virtual bond lengths diff --git a/source/unres/src_MD/parmread.F b/source/unres/src_MD/parmread.F index d6dbc07..68f599a 100644 --- a/source/unres/src_MD/parmread.F +++ b/source/unres/src_MD/parmread.F @@ -211,68 +211,143 @@ C & 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 @@ -281,19 +356,19 @@ 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 @@ -302,8 +377,10 @@ C 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 @@ -651,9 +728,9 @@ C Martix of D parameters for two dimesional fourier series 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) @@ -690,9 +767,9 @@ C Read of Side-chain backbone correlation parameters 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 @@ -703,7 +780,7 @@ cc maxinter is maximum interaction sites 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 @@ -713,7 +790,7 @@ cc maxinter is maximum interaction sites 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 @@ -754,7 +831,7 @@ cc maxinter is maximum interaction sites 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) @@ -768,31 +845,31 @@ cc maxinter is maximum interaction sites 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 @@ -1006,9 +1083,17 @@ C----------------------- 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=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 @@ -1120,7 +1205,7 @@ c augm(i,j)=0.5D0**(2*expon)*aa(i,j) 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 @@ -1158,7 +1243,7 @@ C 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