3 C Read the parameters of the probability distributions of the virtual-bond
4 C valence angles and the side chains and energy parameters.
6 C Important! Energy-term weights ARE NOT read here; they are read from the
7 C main input file instead, because NO defaults have yet been set for these
10 implicit real*8 (a-h,o-z)
16 include 'COMMON.IOUNITS'
17 include 'COMMON.CHAIN'
18 include 'COMMON.INTERACT'
20 include 'COMMON.LOCAL'
21 include 'COMMON.TORSION'
22 include 'COMMON.SCCOR'
23 include 'COMMON.SCROT'
24 include 'COMMON.FFIELD'
25 include 'COMMON.NAMES'
26 include 'COMMON.SBRIDGE'
28 include 'COMMON.SETUP'
30 character*1 onelett(4) /"G","A","P","D"/
32 dimension blower(3,3,maxlob)
34 character*3 lancuch,ucase
36 C For printing parameters after they are read set the following in the UNRES
39 C setenv PRINT_PARM YES
41 C To print parameters in LaTeX format rather than as ASCII tables:
45 call getenv_loc("PRINT_PARM",lancuch)
46 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
47 call getenv_loc("LATEX",lancuch)
48 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
50 dwa16=2.0d0**(1.0d0/6.0d0)
52 C Assign virtual-bond length
57 c Read the virtual-bond parameters, masses, and moments of inertia
58 c and Stokes' radii of the peptide group and side chains
61 read (ibond,*) vbldp0,akp,mp,ip,pstok
64 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
69 dsc_inv(i)=1.0D0/dsc(i)
73 read (ibond,*) junk,vbldp0,akp,rjunk,mp,ip,pstok
75 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
76 & j=1,nbondterm(i)),msc(i),isc(i),restok(i)
81 dsc_inv(i)=1.0D0/dsc(i)
86 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
87 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',
89 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
91 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
92 & vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
94 write (iout,'(13x,3f10.5)')
95 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
101 C Read the parameters of the probability distribution/energy expression
102 C of the virtual-bond valence angles theta
105 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i),j=1,2),
107 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
108 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
109 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
116 & 'Parameters of the virtual-bond valence angles:'
117 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
118 & ' ATHETA0 ',' A1 ',' A2 ',
121 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
122 & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
124 write (iout,'(/a/9x,5a/79(1h-))')
125 & 'Parameters of the expression for sigma(theta_c):',
126 & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
127 & ' ALPH3 ',' SIGMA0C '
129 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
130 & (polthet(j,i),j=0,3),sigc0(i)
132 write (iout,'(/a/9x,5a/79(1h-))')
133 & 'Parameters of the second gaussian:',
134 & ' THETA0 ',' SIGMA0 ',' G1 ',
137 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
138 & sig0(i),(gthet(j,i),j=1,3)
142 & 'Parameters of the virtual-bond valence angles:'
143 write (iout,'(/a/9x,5a/79(1h-))')
144 & 'Coefficients of expansion',
145 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
146 & ' b1*10^1 ',' b2*10^1 '
148 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
149 & a0thet(i),(100*athet(j,i),j=1,2),(10*bthet(j,i),j=1,2)
151 write (iout,'(/a/9x,5a/79(1h-))')
152 & 'Parameters of the expression for sigma(theta_c):',
153 & ' alpha0 ',' alph1 ',' alph2 ',
154 & ' alhp3 ',' sigma0c '
156 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
157 & (polthet(j,i),j=0,3),sigc0(i)
159 write (iout,'(/a/9x,5a/79(1h-))')
160 & 'Parameters of the second gaussian:',
161 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
164 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
165 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
171 C Read the parameters of Utheta determined from ab initio surfaces
172 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
174 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,
175 & ntheterm3,nsingle,ndouble
176 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
177 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
183 aathet(l,i,j,k)=0.0d0
187 bbthet(m,l,i,j,k)=0.0d0
188 ccthet(m,l,i,j,k)=0.0d0
189 ddthet(m,l,i,j,k)=0.0d0
190 eethet(m,l,i,j,k)=0.0d0
196 ffthet(mm,m,l,i,j,k)=0.0d0
197 ggthet(mm,m,l,i,j,k)=0.0d0
207 read (ithep,'(3a)',end=111,err=111) res1,res2,res3
208 read (ithep,*,end=111,err=111) aa0thet(i,j,k)
209 read (ithep,*,end=111,err=111)(aathet(l,i,j,k),l=1,ntheterm)
210 read (ithep,*,end=111,err=111)
211 & ((bbthet(lll,ll,i,j,k),lll=1,nsingle),
212 & (ccthet(lll,ll,i,j,k),lll=1,nsingle),
213 & (ddthet(lll,ll,i,j,k),lll=1,nsingle),
214 & (eethet(lll,ll,i,j,k),lll=1,nsingle),ll=1,ntheterm2)
215 read (ithep,*,end=111,err=111)
216 & (((ffthet(llll,lll,ll,i,j,k),ffthet(lll,llll,ll,i,j,k),
217 & ggthet(llll,lll,ll,i,j,k),ggthet(lll,llll,ll,i,j,k),
218 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
223 C For dummy ends assign glycine-type coefficients of theta-only terms; the
224 C coefficients of theta-and-gamma-dependent terms are zero.
229 aathet(l,i,j,nthetyp+1)=aathet(l,i,j,1)
230 aathet(l,nthetyp+1,i,j)=aathet(l,1,i,j)
232 aa0thet(i,j,nthetyp+1)=aa0thet(i,j,1)
233 aa0thet(nthetyp+1,i,j)=aa0thet(1,i,j)
236 aathet(l,nthetyp+1,i,nthetyp+1)=aathet(l,1,i,1)
238 aa0thet(nthetyp+1,i,nthetyp+1)=aa0thet(1,i,1)
241 C Control printout of the coefficients of virtual-bond-angle potentials
244 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
248 write (iout,'(//4a)')
249 & 'Type ',onelett(i),onelett(j),onelett(k)
250 write (iout,'(//a,10x,a)') " l","a[l]"
251 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k)
252 write (iout,'(i2,1pe15.5)')
253 & (l,aathet(l,i,j,k),l=1,ntheterm)
255 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))')
256 & "b",l,"c",l,"d",l,"e",l
258 write (iout,'(i2,4(1pe15.5))') m,
259 & bbthet(m,l,i,j,k),ccthet(m,l,i,j,k),
260 & ddthet(m,l,i,j,k),eethet(m,l,i,j,k)
264 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
265 & "f+",l,"f-",l,"g+",l,"g-",l
268 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
269 & ffthet(n,m,l,i,j,k),ffthet(m,n,l,i,j,k),
270 & ggthet(n,m,l,i,j,k),ggthet(m,n,l,i,j,k)
279 write (2,*) "Start reading THETA_PDB"
281 read (ithep_pdb,*,err=111,end=111) a0thet(i),(athet(j,i),j=1,2),
283 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
284 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
285 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
288 write (2,*) "End reading THETA_PDB"
294 C Read the parameters of the probability distribution/energy expression
295 C of the side chains.
298 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
302 dsc_inv(i)=1.0D0/dsc(i)
313 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
314 & ((blower(k,l,1),l=1,k),k=1,3)
316 read (irotam,*,end=112,err=112) bsc(j,i)
317 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
318 & ((blower(k,l,j),l=1,k),k=1,3)
325 akl=akl+blower(k,m,j)*blower(l,m,j)
336 write (iout,'(/a)') 'Parameters of side-chain local geometry'
341 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
342 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
343 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
344 & 'log h',(bsc(j,i),j=1,nlobi)
345 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
346 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
348 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
349 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
352 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
353 write (iout,'(a,f10.4,4(16x,f10.4))')
354 & 'Center ',(bsc(j,i),j=1,nlobi)
355 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
364 C Read scrot parameters for potentials determined from all-atom AM1 calculations
365 C added by Urszula Kozlowska 07/11/2007
368 read (irotam,*,end=112,err=112)
370 read (irotam,*,end=112,err=112)
373 read(irotam,*,end=112,err=112) sc_parmin(j,i)
378 C Read the parameters of the probability distribution/energy expression
379 C of the side chains.
382 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
386 dsc_inv(i)=1.0D0/dsc(i)
397 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
398 & ((blower(k,l,1),l=1,k),k=1,3)
400 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
401 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
402 & ((blower(k,l,j),l=1,k),k=1,3)
409 akl=akl+blower(k,m,j)*blower(l,m,j)
424 C Read torsional parameters in old format
426 read (itorp,*,end=113,err=113) ntortyp,nterm_old
427 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
428 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
433 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
439 write (iout,'(/a/)') 'Torsional constants:'
442 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
443 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
449 C Read torsional parameters
451 read (itorp,*,end=113,err=113) ntortyp
452 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
453 c write (iout,*) 'ntortyp',ntortyp
456 read (itorp,*,end=113,err=113) nterm(i,j),nlor(i,j)
460 read (itorp,*,end=113,err=113) kk,v1(k,i,j),v2(k,i,j)
461 v0ij=v0ij+si*v1(k,i,j)
465 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
466 & vlor2(k,i,j),vlor3(k,i,j)
467 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
474 write (iout,'(/a/)') 'Torsional constants:'
477 write (iout,*) 'ityp',i,' jtyp',j
478 write (iout,*) 'Fourier constants'
480 write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
482 write (iout,*) 'Lorenz constants'
484 write (iout,'(3(1pe15.5))')
485 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
491 C 6/23/01 Read parameters for double torsionals
496 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
497 if (t1.ne.onelett(i) .or. t2.ne.onelett(j)
498 & .or. t3.ne.onelett(k)) then
499 write (iout,*) "Error in double torsional parameter file",
502 call MPI_Finalize(Ierror)
504 stop "Error in double torsional parameter file"
506 read (itordp,*,end=114,err=114) ntermd_1(i,j,k),
508 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k),l=1,
510 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k),l=1,
512 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k),l=1,
514 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k),l=1,
516 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k),
517 & v2c(m,l,i,j,k),v2s(l,m,i,j,k),v2s(m,l,i,j,k),
518 & m=1,l-1),l=1,ntermd_2(i,j,k))
524 write (iout,*) 'Constants for double torsionals'
528 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
529 & ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k)
531 write (iout,*) 'Single angles:'
532 do l=1,ntermd_1(i,j,k)
533 write (iout,'(i5,2f10.5,5x,2f10.5)') l,
534 & v1c(1,l,i,j,k),v1s(1,l,i,j,k),
535 & v1c(2,l,i,j,k),v1s(2,l,i,j,k)
538 write (iout,*) 'Pairs of angles:'
539 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
540 do l=1,ntermd_2(i,j,k)
541 write (iout,'(i5,20f10.5)')
542 & l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k))
545 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
546 do l=1,ntermd_2(i,j,k)
547 write (iout,'(i5,20f10.5)')
548 & l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k))
557 C 5/21/07 (AL) Read coefficients of the backbone-local sidechain-local
558 C correlation energies.
560 read (isccor,*,end=119,err=119) nterm_sccor
565 read (isccor,*,end=119,err=119) kk,v1sccor(k,i,j),
572 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
575 write (iout,*) 'ityp',i,' jtyp',j
577 write (iout,'(2(1pe15.5))') v1sccor(k,i,j),v2sccor(k,i,j)
583 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
584 C interaction energy of the Gly, Ala, and Pro prototypes.
588 write (iout,*) "Coefficients of the cumulants"
590 read (ifourier,*) nloctyp
592 read (ifourier,*,end=115,err=115)
593 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
595 write (iout,*) 'Type',i
596 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
622 c Ctilde(1,1,i)=0.0d0
623 c Ctilde(1,2,i)=0.0d0
624 c Ctilde(2,1,i)=0.0d0
625 c Ctilde(2,2,i)=0.0d0
638 c Dtilde(1,1,i)=0.0d0
639 c Dtilde(1,2,i)=0.0d0
640 c Dtilde(2,1,i)=0.0d0
641 c Dtilde(2,2,i)=0.0d0
642 EE(1,1,i)= b(10)+b(11)
643 EE(2,2,i)=-b(10)+b(11)
644 EE(2,1,i)= b(12)-b(13)
645 EE(1,2,i)= b(12)+b(13)
650 c ee(2,1,i)=ee(1,2,i)
654 write (iout,*) 'Type',i
656 write(iout,*) B1(1,i),B1(2,i)
658 write(iout,*) B2(1,i),B2(2,i)
661 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
665 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
669 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
674 C Read electrostatic-interaction parameters
678 write (iout,'(/a)') 'Electrostatic interaction constants:'
679 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
680 & 'IT','JT','APP','BPP','AEL6','AEL3'
682 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
683 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
684 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
685 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
690 app (i,j)=epp(i,j)*rri*rri
691 bpp (i,j)=-2.0D0*epp(i,j)*rri
692 ael6(i,j)=elpp6(i,j)*4.2D0**6
693 ael3(i,j)=elpp3(i,j)*4.2D0**3
694 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
695 & ael6(i,j),ael3(i,j)
699 C Read side-chain interaction parameters.
701 read (isidep,*,end=117,err=117) ipot,expon
702 if (ipot.lt.1 .or. ipot.gt.6) then
703 write (iout,'(2a)') 'Error while reading SC interaction',
704 & 'potential file - unknown potential type.'
706 call MPI_Finalize(Ierror)
712 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
713 & ', exponents are ',expon,2*expon
714 goto (10,20,30,30,40,50) ipot
715 C----------------------- LJ potential ---------------------------------
716 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
717 & (sigma0(i),i=1,ntyp)
719 write (iout,'(/a/)') 'Parameters of the LJ potential:'
720 write (iout,'(a/)') 'The epsilon array:'
721 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
722 write (iout,'(/a)') 'One-body parameters:'
723 write (iout,'(a,4x,a)') 'residue','sigma'
724 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
727 C----------------------- LJK potential --------------------------------
728 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
729 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
731 write (iout,'(/a/)') 'Parameters of the LJK potential:'
732 write (iout,'(a/)') 'The epsilon array:'
733 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
734 write (iout,'(/a)') 'One-body parameters:'
735 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
736 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
740 C---------------------- GB or BP potential -----------------------------
741 30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
742 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
744 C For the GB potential convert sigma'**2 into chi'
747 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
751 write (iout,'(/a/)') 'Parameters of the BP potential:'
752 write (iout,'(a/)') 'The epsilon array:'
753 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
754 write (iout,'(/a)') 'One-body parameters:'
755 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
757 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
758 & chip(i),alp(i),i=1,ntyp)
761 C--------------------- GBV potential -----------------------------------
762 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
763 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
764 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
766 write (iout,'(/a/)') 'Parameters of the GBV potential:'
767 write (iout,'(a/)') 'The epsilon array:'
768 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
769 write (iout,'(/a)') 'One-body parameters:'
770 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
771 & 's||/s_|_^2',' chip ',' alph '
772 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
773 & sigii(i),chip(i),alp(i),i=1,ntyp)
777 C--------------------- Momo potential -----------------------------------
781 read (isidep,*) (icharge(i),i=1,ntyp)
782 c write (2,*) "icharge",(icharge(i),i=1,ntyp)
785 c! write (*,*) "Im in ", i, " ", j
787 & eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i),
788 & (alphasur(k,i,j),k=1,4),sigmap(i,j),sigmap(j,i),
789 & chis(i,j),chis(j,i),
790 & nstate(i,j),(wstate(k,i,j),k=1,4),
791 & ((dhead(l,k,i,j),k=1,2),l=1,2),dtail(i,j),dtail(j,i),
792 & epshead(i,j),sig0head(i,j),
793 & rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j),
794 & alphapol(i,j),alphapol(j,i),
795 & (alphiso(k,i,j),k=1,4),sigiso(i,j),sigiso(j,i),epsintab(i,j)
799 c! write (*,*) "Parameters read in"
801 c! THIS LOOP FILLS PARAMETERS FOR PAIRS OF AA's NOT EXPLICITLY
802 c! DEFINED IN SCPARM.MOMO. IT DOES SO BY TAKING THEM FROM SYMMETRIC
803 c! PAIR, FOR EG. IF ARG-HIS IS BLANK, IT WILL TRY TO TAKE PARAMETERS
806 c! DISABLE IT AT >>YOUR OWN PERIL<<
811 sigma(i,j) = sigma(j,i)
812 nstate(i,j) = nstate(j,i)
815 alphasur(k,i,j) = alphasur(k,j,i)
816 wstate(k,i,j) = wstate(k,j,i)
817 alphiso(k,i,j) = alphiso(k,j,i)
822 dhead(k,l,i,j) = dhead(k,l,j,i)
826 epshead(i,j) = epshead(j,i)
827 sig0head(i,j) = sig0head(j,i)
830 wqdip(k,i,j) = wqdip(k,j,i)
833 wquad(i,j) = wquad(j,i)
834 epsintab(i,j) = epsintab(j,i)
842 c! write (*,'(f17.15)') ((eps(i,j),j=i,ntyp), i=1, ntyp)
843 c! write (*,'(f17.15)') (sigma0(i),i=1,ntyp)
844 c! write (*,'(f17.15)') (sigii(i),i=1,ntyp)
845 c! write (*,'(f17.15)') (chip(i),i=1,ntyp)
846 c! write (*,'(f17.15)') (alp(i),i=1,ntyp)
851 c! This is older parameter-filling loop
855 c IF ((sig0head(i,j).EQ.0.0d0).AND.(sig0head(j,i).NE.0.0d0)) THEN
856 c sig0head(i,j) = sig0head(j,i)
857 c ELSE IF ((sig0head(j,i).EQ.0.0d0).AND.(sig0head(i,j).NE.0.0d0))
859 c sig0head(j,i) = sig0head(i,j)
862 c IF ((epshead(i,j).EQ.0.0d0).AND.(epshead(j,i).NE.0.0d0)) THEN
863 c epshead(i,j) = epshead(j,i)
864 c ELSE IF ((epshead(j,i).EQ.0.0d0).AND.(epshead(i,j).NE.0.0d0))
866 c epshead(j,i) = epshead(i,j)
869 c IF ((wquad(i,j).EQ.0.0d0).AND.(wquad(j,i).NE.0.0d0)) THEN
870 c wquad(i,j) = wquad(j,i)
871 c ELSE IF ((wquad(j,i).EQ.0.0d0).AND.(wquad(i,j).NE.0.0d0))
873 c wquad(j,i) = wquad(i,j)
876 c IF ((sigiso(i,j).EQ.0.0d0).AND.sigiso(j,i).NE.0.0d0)) THEN
877 c sigiiso(i,j) = sigiso(j,i)
878 c ELSE IF ((epshead(j,i).EQ.0.0d0).AND.(epshead(i,j).NE.0.0d0))
880 c sigiso(j,i) = sigiso(i,j)
883 c! IF ((wstate(k,i,j).EQ.0.0d0).AND.
884 c! & (wstate(k,j,i).NE.0.0d0)) THEN
885 c! wstate(k,i,j) = wstate(k,j,i)
886 c! ELSE IF ((wstate(k,j,i).EQ.0.0d0).AND.
887 c! & (wstate(k,i,j).NE.0.0d0)) THEN
888 c! wstate(k,j,i) = wstate(k,i,j)
893 c! THE LOOP MAKING MOMO_TABLES SYMMETRIC ENDS HERE
895 if (.not.lprint) goto 70
897 & "Parameters of the new physics-based SC-SC interaction potential"
898 write (iout,'(/7a)') 'Residues',' epsGB',' rGB',
899 & ' chi1GB',' chi2GB',' chip1GB',' chip2GB'
902 write (iout,'(2(a3,1x),1pe10.3,5(0pf10.3))')
903 & restyp(i),restyp(j),eps(i,j),sigma(i,j),chi(i,j),chi(j,i),
904 & chipp(i,j),chipp(j,i)
907 write (iout,'(/9a)') 'Residues',' alphasur1',' alphasur2',
908 & ' alphasur3',' alphasur4',' sigmap1',' sigmap2',
912 write (iout,'(2(a3,1x),8(0pf10.3))')
913 & restyp(i),restyp(j),(alphasur(k,i,j),k=1,4),
914 & sigmap(i,j),sigmap(j,i),chis(i,j),chis(j,i)
917 write (iout,'(/14a)') 'Residues',' nst ',' wst1',
918 & ' wst2',' wst3',' wst4',' dh11',' dh21',
919 & ' dh12',' dh22',' dt1',' dt2',' epsh1',
923 write (iout,'(2(a3,1x),i3,4f8.3,6f7.2,f9.5,f7.2)')
924 & restyp(i),restyp(j),nstate(i,j),(wstate(k,i,j),k=1,4),
925 & ((dhead(l,k,i,j),l=1,2),k=1,2),dtail(i,j),dtail(j,i),
926 & epshead(i,j),sig0head(i,j)
929 write (iout,'(/12a)') 'Residues',' ch1',' ch2',
930 & ' rborn1',' rborn2',' wqdip1',' wqdip2',
934 write (iout,'(2(a3,1x),2i4,5f10.3)')
935 & restyp(i),restyp(j),icharge(i),icharge(j),
936 & rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j)
939 write (iout,'(/12a)') 'Residues',
941 & ' alphpol2',' alphiso1',' alpiso2',
942 & ' alpiso3',' alpiso4',' sigiso1',' sigiso2',
946 write (iout,'(2(a3,1x),11f10.3)')
947 & restyp(i),restyp(j),alphapol(i,j),alphapol(j,i),
948 & (alphiso(k,i,j),k=1,4),sigiso(i,j),sigiso(j,i),epsintab(i,j)
953 C For the GB potential convert sigma'**2 into chi'
955 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
958 write (iout,'(/a/)') 'Parameters of the BP potential:'
959 write (iout,'(a/)') 'The epsilon array:'
960 CALL printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
961 write (iout,'(/a)') 'One-body parameters:'
962 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
964 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
965 & chip(i),alp(i),i=1,ntyp)
970 C-----------------------------------------------------------------------
971 C Calculate the "working" parameters of SC interactions.
976 alphasur(k,j,i)=alphasur(k,i,j)
977 alphiso(k,j,i)=alphiso(k,i,j)
978 wstate(k,j,i)=wstate(k,i,j)
981 wqdip(k,j,i)=wqdip(k,i,j)
985 dhead(l,k,j,i)=dhead(l,k,i,j)
993 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
994 sigma(j,i)=sigma(i,j)
995 rs0(i,j)=dwa16*sigma(i,j)
1003 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1004 & 'Working parameters of the SC interactions:',
1005 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1010 IF (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4 .or. ipot.eq.6)
1020 sigeps=dsign(1.0D0,epsij)
1022 aa(i,j)=epsij*rrij*rrij
1023 bb(i,j)=-sigeps*epsij*rrij
1026 IF ((ipot.GT.2).AND.(ipot.LT.6)) THEN
1027 sigt1sq = sigma0(i)**2
1028 sigt2sq = sigma0(j)**2
1031 ratsig1 = sigt2sq/sigt1sq
1032 ratsig2 = 1.0D0/ratsig1
1033 chi(i,j) = (sigii1-1.0D0) / (sigii1+ratsig1)
1034 IF (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1035 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1039 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1040 sigmaii(i,j)=rsum_max
1041 sigmaii(j,i)=rsum_max
1043 c sigmaii(i,j)=r0(i,j)
1044 c sigmaii(j,i)=r0(i,j)
1046 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1047 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1048 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1049 augm(i,j)=epsij*r_augm**(2*expon)
1050 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1059 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1060 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1061 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1063 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3),2i3,10f8.4,
1065 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1066 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i),
1067 & icharge(i),icharge(j),chipp(i,j),chipp(j,i),
1068 & (alphasur(k,i,j),k=1,4),sigmap(i,j),sigmap(j,i),
1069 & chis(i,j),chis(j,i),
1070 & nstate(i,j),(wstate(k,i,j),k=1,4),
1071 & ((dhead(l,k,i,j),l=1,2),k=1,2),dtail(i,j),dtail(j,i),
1072 & epshead(i,j),sig0head(i,j),
1073 & rborn(i,j),(wqdip(k,i,j),k=1,2),wquad(i,j),
1074 & alphapol(i,j),alphapol(j,i),
1075 & (alphiso(k,i,j),k=1,4),sigiso(i,j)
1086 C Define the SC-p interaction constants (hard-coded; old style)
1089 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1091 c aad(i,1)=0.3D0*4.0D0**12
1092 C Following line for constants currently implemented
1093 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1094 aad(i,1)=1.5D0*4.0D0**12
1095 c aad(i,1)=0.17D0*5.6D0**12
1097 C "Soft" SC-p repulsion
1099 C Following line for constants currently implemented
1100 c aad(i,1)=0.3D0*4.0D0**6
1101 C "Hard" SC-p repulsion
1102 bad(i,1)=3.0D0*4.0D0**6
1103 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1112 C 8/9/01 Read the SC-p interaction constants from file
1115 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1118 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1119 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1120 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1121 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1125 write (iout,'(/a)') "Parameters of SC-p interactions:"
1127 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1128 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1133 C Define the constants of the disulfide bridge
1137 c Old arbitrary potential - commented out.
1142 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1143 c energy surface of diethyl disulfide.
1144 c A. Liwo and U. Kozlowska, 11/24/03
1161 c! print *, "END OF PARMREAD"
1162 c! print *, "eps = ", eps(9,9), "sigma = ", sigma(9,9)
1163 c! print *, "chi1 = ", chi(9,9), "chi2 = ", chi(9,9)
1164 c! print *, "chip1 = ", chipp(9,9), "chip2 = ", chipp(9,9)
1165 c! print *, "sig1 = ", sigmap(9,9), "sig2 = ", sigmap(9,9)
1166 c! print *, "chis1 = ", chis(9,9)," chis2 = ", chis(9,9)
1167 c! print *, "END OF PARMREAD"
1172 write (iout,'(/a)') "Disulfide bridge parameters:"
1173 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1174 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1175 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1176 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1180 111 write (iout,*) "Error reading bending energy parameters."
1182 112 write (iout,*) "Error reading rotamer energy parameters."
1184 113 write (iout,*) "Error reading torsional energy parameters."
1186 114 write (iout,*) "Error reading double torsional energy parameters."
1189 & "Error reading cumulant (multibody energy) parameters."
1191 116 write (iout,*) "Error reading electrostatic energy parameters."
1193 117 write (iout,*) "Error reading side chain interaction parameters."
1195 118 write (iout,*) "Error reading SCp interaction parameters."
1197 119 write (iout,*) "Error reading SCCOR parameters"
1200 call MPI_Finalize(Ierror)
1207 subroutine getenv_loc(var, val)
1208 character(*) var, val
1211 character(2000) line
1214 open (196,file='env',status='old',readonly,shared)
1216 c write(*,*)'looking for ',var
1217 10 read(196,*,err=11,end=11)line
1218 iread=index(line,var)
1219 c write(*,*)iread,' ',var,' ',line
1220 if (iread.eq.0) go to 10
1221 c write(*,*)'---> ',line
1227 iread=iread+ilen(var)+1
1228 read (line(iread:),*,err=12,end=12) val
1229 c write(*,*)'OK: ',var,' = ',val
1235 #elif (defined CRAY)
1236 integer lennam,lenval,ierror
1238 c getenv using a POSIX call, useful on the T3D
1239 c Sept 1996, comment out error check on advice of H. Pritchard
1242 if(lennam.le.0) stop '--error calling getenv--'
1243 call pxfgetenv(var,lennam,val,lenval,ierror)
1244 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1246 call getenv(var,val)