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))
556 C Read of Side-chain backbone correlation parameters
557 C Modified 11 May 2012 by Adasko
560 read (isccor,*,end=1113,err=1113) nsccortyp
561 read (isccor,*,end=1113,err=1113) (isccortyp(i),i=1,ntyp)
562 write (iout,*) 'ntortyp',ntortyp
564 cc maxinter is maximum interaction sites
568 read (isccor,*,end=1113,err=1113) nterm_sccor(i,j),
573 DO k = 1, nterm_sccor(i,j)
574 read (isccor,*,end=1113,err=1113) kk, v1sccor(k,l,i,j)
576 v0ijsccor = v0ijsccor + si * v1sccor(k,l,i,j)
579 DO k = 1, nlor_sccor(i,j)
580 read (isccor,*,end=1113,err=1113) kk, vlor1sccor(k,i,j),
581 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
582 v0ijsccor = v0ijsccor + vlor1sccor(k,i,j) /
583 & (1 + vlor3sccor(k,i,j)**2)
585 v0sccor(i,j) = v0ijsccor
592 write (iout,'(/a/)') 'Torsional constants:'
595 write (iout,*) 'ityp',i,' jtyp',j
596 write (iout,*) 'Fourier constants'
597 DO k=1,nterm_sccor(i,j)
598 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
600 write (iout,*) 'Lorenz constants'
601 DO k=1,nlor_sccor(i,j)
602 write (iout,'(3(1pe15.5))')
603 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
610 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
611 C interaction energy of the Gly, Ala, and Pro prototypes.
615 write (iout,*) "Coefficients of the cumulants"
617 read (ifourier,*) nloctyp
619 read (ifourier,*,end=115,err=115)
620 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
622 write (iout,*) 'Type',i
623 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
649 c Ctilde(1,1,i)=0.0d0
650 c Ctilde(1,2,i)=0.0d0
651 c Ctilde(2,1,i)=0.0d0
652 c Ctilde(2,2,i)=0.0d0
665 c Dtilde(1,1,i)=0.0d0
666 c Dtilde(1,2,i)=0.0d0
667 c Dtilde(2,1,i)=0.0d0
668 c Dtilde(2,2,i)=0.0d0
669 EE(1,1,i)= b(10)+b(11)
670 EE(2,2,i)=-b(10)+b(11)
671 EE(2,1,i)= b(12)-b(13)
672 EE(1,2,i)= b(12)+b(13)
677 c ee(2,1,i)=ee(1,2,i)
681 write (iout,*) 'Type',i
683 write(iout,*) B1(1,i),B1(2,i)
685 write(iout,*) B2(1,i),B2(2,i)
688 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
692 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
696 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
701 C Read electrostatic-interaction parameters
705 write (iout,'(/a)') 'Electrostatic interaction constants:'
706 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
707 & 'IT','JT','APP','BPP','AEL6','AEL3'
709 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
710 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
711 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
712 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
717 app (i,j)=epp(i,j)*rri*rri
718 bpp (i,j)=-2.0D0*epp(i,j)*rri
719 ael6(i,j)=elpp6(i,j)*4.2D0**6
720 ael3(i,j)=elpp3(i,j)*4.2D0**3
721 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
722 & ael6(i,j),ael3(i,j)
726 C Read side-chain interaction parameters.
728 read (isidep,*,end=117,err=117) ipot,expon
729 if (ipot.lt.1 .or. ipot.gt.6) then
730 write (iout,'(2a)') 'Error while reading SC interaction',
731 & 'potential file - unknown potential type.'
733 call MPI_Finalize(Ierror)
739 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
740 & ', exponents are ',expon,2*expon
741 goto (10,20,30,30,40,50) ipot
742 C----------------------- LJ potential ---------------------------------
743 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
744 & (sigma0(i),i=1,ntyp)
746 write (iout,'(/a/)') 'Parameters of the LJ potential:'
747 write (iout,'(a/)') 'The epsilon array:'
748 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
749 write (iout,'(/a)') 'One-body parameters:'
750 write (iout,'(a,4x,a)') 'residue','sigma'
751 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
754 C----------------------- LJK potential --------------------------------
755 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
756 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
758 write (iout,'(/a/)') 'Parameters of the LJK potential:'
759 write (iout,'(a/)') 'The epsilon array:'
760 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
761 write (iout,'(/a)') 'One-body parameters:'
762 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
763 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
767 C---------------------- GB or BP potential -----------------------------
768 30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
769 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
771 C For the GB potential convert sigma'**2 into chi'
774 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
778 write (iout,'(/a/)') 'Parameters of the GB or BP potential:'
779 write (iout,'(a/)') 'The epsilon array:'
780 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
781 write (iout,'(/a)') 'One-body parameters:'
782 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
784 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
785 & chip(i),alp(i),i=1,ntyp)
788 C--------------------- GBV potential -----------------------------------
789 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
790 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
791 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
793 write (iout,'(/a/)') 'Parameters of the GBV potential:'
794 write (iout,'(a/)') 'The epsilon array:'
795 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
796 write (iout,'(/a)') 'One-body parameters:'
797 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
798 & 's||/s_|_^2',' chip ',' alph '
799 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
800 & sigii(i),chip(i),alp(i),i=1,ntyp)
804 C--------------------- Momo potential -----------------------------------
808 read (isidep,*) (icharge(i),i=1,ntyp)
809 c write (2,*) "icharge",(icharge(i),i=1,ntyp)
812 c! write (*,*) "Im in ", i, " ", j
814 & eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i),
815 & (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j),
816 & chis(i,j),chis(j,i),
817 & nstate(i,j),(wstate(k,i,j),k=1,4),
822 & dtail(1,i,j),dtail(2,i,j),
823 & epshead(i,j),sig0head(i,j),
824 & rborn(i,j),rborn(j,i),
825 & (wqdip(k,i,j),k=1,2),wquad(i,j),
826 & alphapol(i,j),alphapol(j,i),
827 & (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j)
830 c! write (*,*) "nstate gly-gly", nstate(10,10)
831 c! THIS LOOP FILLS PARAMETERS FOR PAIRS OF AA's NOT EXPLICITLY
832 c! DEFINED IN SCPARM.MOMO. IT DOES SO BY TAKING THEM FROM SYMMETRIC
833 c! PAIR, FOR EG. IF ARG-HIS IS BLANK, IT WILL TAKE PARAMETERS
836 c! DISABLE IT AT >>YOUR OWN PERIL<<
841 sigma(i,j) = sigma(j,i)
842 nstate(i,j) = nstate(j,i)
843 sigmap1(i,j) = sigmap1(j,i)
844 sigmap2(i,j) = sigmap2(j,i)
845 sigiso1(i,j) = sigiso1(j,i)
846 sigiso2(i,j) = sigiso2(j,i)
849 alphasur(k,i,j) = alphasur(k,j,i)
850 wstate(k,i,j) = wstate(k,j,i)
851 alphiso(k,i,j) = alphiso(k,j,i)
854 dhead(2,1,i,j) = dhead(1,1,j,i)
855 dhead(2,2,i,j) = dhead(1,2,j,i)
856 dhead(1,1,i,j) = dhead(2,1,j,i)
857 dhead(1,2,i,j) = dhead(2,2,j,i)
858 dtail(1,i,j) = dtail(1,j,i)
859 dtail(2,i,j) = dtail(2,j,i)
862 c! write (*,*) "dhead(k,l,j,i) = ", dhead(k,l,j,i)
863 c! write (*,*) "dhead(k,l,i,j) = ", dhead(k,l,i,j)
864 c! dhead(l,k,i,j) = dhead(k,l,j,i)
868 epshead(i,j) = epshead(j,i)
869 sig0head(i,j) = sig0head(j,i)
872 wqdip(k,i,j) = wqdip(k,j,i)
875 wquad(i,j) = wquad(j,i)
876 epsintab(i,j) = epsintab(j,i)
881 if (.not.lprint) goto 70
883 & "Parameters of the new physics-based SC-SC interaction potential"
884 write (iout,'(/7a)') 'Residues',' epsGB',' rGB',
885 & ' chi1GB',' chi2GB',' chip1GB',' chip2GB'
888 write (iout,'(2(a3,1x),1pe10.3,5(0pf10.3))')
889 & restyp(i),restyp(j),eps(i,j),sigma(i,j),chi(i,j),chi(j,i),
890 & chipp(i,j),chipp(j,i)
893 write (iout,'(/9a)') 'Residues',' alphasur1',' alphasur2',
894 & ' alphasur3',' alphasur4',' sigmap1',' sigmap2',
898 write (iout,'(2(a3,1x),8(0pf10.3))')
899 & restyp(i),restyp(j),(alphasur(k,i,j),k=1,4),
900 & sigmap1(i,j),sigmap2(j,i),chis(i,j),chis(j,i)
903 write (iout,'(/14a)') 'Residues',' nst ',' wst1',
904 & ' wst2',' wst3',' wst4',' dh11',' dh21',
905 & ' dh12',' dh22',' dt1',' dt2',' epsh1',
909 write (iout,'(2(a3,1x),i3,4f8.3,6f7.2,f9.5,f7.2)')
910 & restyp(i),restyp(j),nstate(i,j),(wstate(k,i,j),k=1,4),
911 & ((dhead(l,k,i,j),l=1,2),k=1,2),dtail(1,i,j),dtail(2,i,j),
912 & epshead(i,j),sig0head(i,j)
915 write (iout,'(/12a)') 'Residues',' ch1',' ch2',
916 & ' rborn1',' rborn2',' wqdip1',' wqdip2',
920 write (iout,'(2(a3,1x),2i4,5f10.3)')
921 & restyp(i),restyp(j),icharge(i),icharge(j),
922 & rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j)
925 write (iout,'(/12a)') 'Residues',
927 & ' alphpol2',' alphiso1',' alpiso2',
928 & ' alpiso3',' alpiso4',' sigiso1',' sigiso2',
932 write (iout,'(2(a3,1x),11f10.3)')
933 & restyp(i),restyp(j),alphapol(i,j),alphapol(j,i),
934 & (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(j,i),
940 C For the GB potential convert sigma'**2 into chi'
942 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
945 write (iout,'(/a/)') 'Parameters of the BP potential:'
946 write (iout,'(a/)') 'The epsilon array:'
947 CALL printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
948 write (iout,'(/a)') 'One-body parameters:'
949 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
951 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
952 & chip(i),alp(i),i=1,ntyp)
957 C-----------------------------------------------------------------------
958 C Calculate the "working" parameters of SC interactions.
963 alphasur(k,j,i)=alphasur(k,i,j)
964 alphiso(k,j,i)=alphiso(k,i,j)
965 wstate(k,j,i)=wstate(k,i,j)
968 wqdip(k,j,i)=wqdip(k,i,j)
972 dhead(l,k,j,i)=dhead(l,k,i,j)
980 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
981 sigma(j,i)=sigma(i,j)
982 rs0(i,j)=dwa16*sigma(i,j)
989 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
990 & 'Working parameters of the SC interactions:',
991 & ' a ',' b ',' augm ',' sigma ',' r0 ',
996 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4 .or. ipot.eq.6 ) THEN
1005 sigeps=dsign(1.0D0,epsij)
1007 aa(i,j)=epsij*rrij*rrij
1008 bb(i,j)=-sigeps*epsij*rrij
1011 IF ((ipot.gt.2).AND.(ipot.LT.6)) THEN
1012 sigt1sq=sigma0(i)**2
1013 sigt2sq=sigma0(j)**2
1016 ratsig1=sigt2sq/sigt1sq
1017 ratsig2=1.0D0/ratsig1
1018 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1019 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1020 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1024 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1025 sigmaii(i,j)=rsum_max
1026 sigmaii(j,i)=rsum_max
1028 c sigmaii(i,j)=r0(i,j)
1029 c sigmaii(j,i)=r0(i,j)
1031 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1032 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1033 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1034 augm(i,j)=epsij*r_augm**(2*expon)
1035 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1043 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1044 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1045 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1047 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3),2i3,10f8.4,
1049 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1050 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i),
1051 & icharge(i),icharge(j),chipp(i,j),chipp(j,i),
1052 & (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(j,i),
1053 & chis(i,j),chis(j,i),
1054 & nstate(i,j),(wstate(k,i,j),k=1,4),
1055 & ((dhead(l,k,i,j),l=1,2),k=1,2),dtail(1,i,j),dtail(2,i,j),
1056 & epshead(i,j),sig0head(i,j),
1057 & rborn(i,j),(wqdip(k,i,j),k=1,2),wquad(i,j),
1058 & alphapol(i,j),alphapol(j,i),
1059 & (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j)
1067 C Define the SC-p interaction constants (hard-coded; old style)
1070 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1072 c aad(i,1)=0.3D0*4.0D0**12
1073 C Following line for constants currently implemented
1074 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1075 aad(i,1)=1.5D0*4.0D0**12
1076 c aad(i,1)=0.17D0*5.6D0**12
1078 C "Soft" SC-p repulsion
1080 C Following line for constants currently implemented
1081 c aad(i,1)=0.3D0*4.0D0**6
1082 C "Hard" SC-p repulsion
1083 bad(i,1)=3.0D0*4.0D0**6
1084 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1093 C 8/9/01 Read the SC-p interaction constants from file
1096 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1099 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1100 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1101 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1102 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1106 write (iout,*) "Parameters of SC-p interactions:"
1108 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1109 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1114 C Define the constants of the disulfide bridge
1118 c Old arbitrary potential - commented out.
1123 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1124 c energy surface of diethyl disulfide.
1125 c A. Liwo and U. Kozlowska, 11/24/03
1142 write (iout,'(/a)') "Disulfide bridge parameters:"
1143 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1144 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1145 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1146 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1150 111 write (iout,*) "Error reading bending energy parameters."
1152 112 write (iout,*) "Error reading rotamer energy parameters."
1154 113 write (iout,*) "Error reading torsional energy parameters."
1157 & "Error reading side-chain torsional energy parameters."
1159 114 write (iout,*) "Error reading double torsional energy parameters."
1162 & "Error reading cumulant (multibody energy) parameters."
1164 116 write (iout,*) "Error reading electrostatic energy parameters."
1166 117 write (iout,*) "Error reading side chain interaction parameters."
1168 118 write (iout,*) "Error reading SCp interaction parameters."
1170 119 write (iout,*) "Error reading SCCOR parameters"
1173 call MPI_Finalize(Ierror)
1180 subroutine getenv_loc(var, val)
1181 character(*) var, val
1184 character(2000) line
1187 open (196,file='env',status='old',readonly,shared)
1189 c write(*,*)'looking for ',var
1190 10 read(196,*,err=11,end=11)line
1191 iread=index(line,var)
1192 c write(*,*)iread,' ',var,' ',line
1193 if (iread.eq.0) go to 10
1194 c write(*,*)'---> ',line
1200 iread=iread+ilen(var)+1
1201 read (line(iread:),*,err=12,end=12) val
1202 c write(*,*)'OK: ',var,' = ',val
1208 #elif (defined CRAY)
1209 integer lennam,lenval,ierror
1211 c getenv using a POSIX call, useful on the T3D
1212 c Sept 1996, comment out error check on advice of H. Pritchard
1215 if(lennam.le.0) stop '--error calling getenv--'
1216 call pxfgetenv(var,lennam,val,lenval,ierror)
1217 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1219 call getenv(var,val)