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.
586 read (ifourier,*,end=115,err=115) nloctyp
589 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
591 read (ifourier,*,end=115,err=115) (bnew1(ii,1,i),ii=1,3)
592 read (ifourier,*,end=115,err=115) (bnew2(ii,1,i),ii=1,3)
593 read (ifourier,*,end=115,err=115) (bnew1(ii,2,i),ii=1,1)
594 read (ifourier,*,end=115,err=115) (bnew2(ii,2,i),ii=1,1)
595 read (ifourier,*,end=115,err=115) (eenew(ii,i),ii=1,1)
598 write (iout,*) 'Type',i
599 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
601 write (iout,'("Bnew1",i3,2f10.5)') ii,(bnew1(ii,k,i),k=1,2)
602 write (iout,'("Bnew2",i3,2f10.5)') ii,(bnew2(ii,k,i),k=1,2)
605 write (iout,'("EEnew",i3,2f10.5)') ii,eenew(ii,i)
611 B1tilde(1,i) = b(3,i)
612 B1tilde(2,i) =-b(5,i)
622 Ctilde(2,1,i)=-b(9,i)
630 Dtilde(2,1,i)=-b(8,i)
633 EEold(1,1,i)= b(10,i)+b(11,i)
634 EEold(2,2,i)=-b(10,i)+b(11,i)
635 EEold(2,1,i)= b(12,i)-b(13,i)
636 EEold(1,2,i)= b(12,i)+b(13,i)
637 c EEold(1,1,-i)= b(10,i)+b(11,i)
638 c EEold(2,2,-i)=-b(10,i)+b(11,i)
639 c EEold(2,1,-i)=-b(12,i)+b(13,i)
640 c EEold(1,2,-i)=-b(12,i)-b(13,i)
642 EE(1,1,i)= b(10,i)+b(11,i)
643 EE(2,2,i)=-b(10,i)+b(11,i)
644 EE(2,1,i)= b(12,i)-b(13,i)
645 EE(1,2,i)= b(12,i)+b(13,i)
649 write (iout,*) 'Type',i
651 write(iout,*) B1(1,i),B1(2,i)
653 write(iout,*) B2(1,i),B2(2,i)
656 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
660 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
664 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
669 C Read electrostatic-interaction parameters
673 write (iout,'(/a)') 'Electrostatic interaction constants:'
674 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
675 & 'IT','JT','APP','BPP','AEL6','AEL3'
677 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
678 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
679 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
680 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
685 app (i,j)=epp(i,j)*rri*rri
686 bpp (i,j)=-2.0D0*epp(i,j)*rri
687 ael6(i,j)=elpp6(i,j)*4.2D0**6
688 ael3(i,j)=elpp3(i,j)*4.2D0**3
689 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
690 & ael6(i,j),ael3(i,j)
694 C Read side-chain interaction parameters.
696 read (isidep,*,end=117,err=117) ipot,expon
697 if (ipot.lt.1 .or. ipot.gt.6) then
698 write (iout,'(2a)') 'Error while reading SC interaction',
699 & 'potential file - unknown potential type.'
701 call MPI_Finalize(Ierror)
707 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
708 & ', exponents are ',expon,2*expon
709 goto (10,20,30,30,40,50) ipot
710 C----------------------- LJ potential ---------------------------------
711 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
712 & (sigma0(i),i=1,ntyp)
714 write (iout,'(/a/)') 'Parameters of the LJ potential:'
715 write (iout,'(a/)') 'The epsilon array:'
716 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
717 write (iout,'(/a)') 'One-body parameters:'
718 write (iout,'(a,4x,a)') 'residue','sigma'
719 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
722 C----------------------- LJK potential --------------------------------
723 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
724 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
726 write (iout,'(/a/)') 'Parameters of the LJK potential:'
727 write (iout,'(a/)') 'The epsilon array:'
728 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
729 write (iout,'(/a)') 'One-body parameters:'
730 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
731 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
735 C---------------------- GB or BP potential -----------------------------
736 30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
737 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
739 C For the GB potential convert sigma'**2 into chi'
742 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
746 write (iout,'(/a/)') 'Parameters of the BP 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,4a)') 'residue',' sigma ','s||/s_|_^2',
752 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
753 & chip(i),alp(i),i=1,ntyp)
756 C--------------------- GBV potential -----------------------------------
757 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
758 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
759 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
761 write (iout,'(/a/)') 'Parameters of the GBV potential:'
762 write (iout,'(a/)') 'The epsilon array:'
763 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
764 write (iout,'(/a)') 'One-body parameters:'
765 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
766 & 's||/s_|_^2',' chip ',' alph '
767 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
768 & sigii(i),chip(i),alp(i),i=1,ntyp)
772 C--------------------- Momo potential -----------------------------------
776 read (isidep,*) (icharge(i),i=1,ntyp)
777 c write (2,*) "icharge",(icharge(i),i=1,ntyp)
780 c! write (*,*) "Im in ", i, " ", j
782 & eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i),
783 & (alphasur(k,i,j),k=1,4),sigmap(i,j),sigmap(j,i),
784 & chis(i,j),chis(j,i),
785 & nstate(i,j),(wstate(k,i,j),k=1,4),
786 & ((dhead(l,k,i,j),k=1,2),l=1,2),dtail(i,j),dtail(j,i),
787 & epshead(i,j),sig0head(i,j),
788 & rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j),
789 & alphapol(i,j),alphapol(j,i),
790 & (alphiso(k,i,j),k=1,4),sigiso(i,j),sigiso(j,i),epsintab(i,j)
794 c! write (*,*) "Parameters read in"
796 c! THIS LOOP FILLS PARAMETERS FOR PAIRS OF AA's NOT EXPLICITLY
797 c! DEFINED IN SCPARM.MOMO. IT DOES SO BY TAKING THEM FROM SYMMETRIC
798 c! PAIR, FOR EG. IF ARG-HIS IS BLANK, IT WILL TRY TO TAKE PARAMETERS
801 c! DISABLE IT AT >>YOUR OWN PERIL<<
806 sigma(i,j) = sigma(j,i)
807 nstate(i,j) = nstate(j,i)
810 alphasur(k,i,j) = alphasur(k,j,i)
811 wstate(k,i,j) = wstate(k,j,i)
812 alphiso(k,i,j) = alphiso(k,j,i)
817 dhead(k,l,i,j) = dhead(k,l,j,i)
821 epshead(i,j) = epshead(j,i)
822 sig0head(i,j) = sig0head(j,i)
825 wqdip(k,i,j) = wqdip(k,j,i)
828 wquad(i,j) = wquad(j,i)
829 epsintab(i,j) = epsintab(j,i)
837 c! write (*,'(f17.15)') ((eps(i,j),j=i,ntyp), i=1, ntyp)
838 c! write (*,'(f17.15)') (sigma0(i),i=1,ntyp)
839 c! write (*,'(f17.15)') (sigii(i),i=1,ntyp)
840 c! write (*,'(f17.15)') (chip(i),i=1,ntyp)
841 c! write (*,'(f17.15)') (alp(i),i=1,ntyp)
846 c! This is older parameter-filling loop
850 c IF ((sig0head(i,j).EQ.0.0d0).AND.(sig0head(j,i).NE.0.0d0)) THEN
851 c sig0head(i,j) = sig0head(j,i)
852 c ELSE IF ((sig0head(j,i).EQ.0.0d0).AND.(sig0head(i,j).NE.0.0d0))
854 c sig0head(j,i) = sig0head(i,j)
857 c IF ((epshead(i,j).EQ.0.0d0).AND.(epshead(j,i).NE.0.0d0)) THEN
858 c epshead(i,j) = epshead(j,i)
859 c ELSE IF ((epshead(j,i).EQ.0.0d0).AND.(epshead(i,j).NE.0.0d0))
861 c epshead(j,i) = epshead(i,j)
864 c IF ((wquad(i,j).EQ.0.0d0).AND.(wquad(j,i).NE.0.0d0)) THEN
865 c wquad(i,j) = wquad(j,i)
866 c ELSE IF ((wquad(j,i).EQ.0.0d0).AND.(wquad(i,j).NE.0.0d0))
868 c wquad(j,i) = wquad(i,j)
871 c IF ((sigiso(i,j).EQ.0.0d0).AND.sigiso(j,i).NE.0.0d0)) THEN
872 c sigiiso(i,j) = sigiso(j,i)
873 c ELSE IF ((epshead(j,i).EQ.0.0d0).AND.(epshead(i,j).NE.0.0d0))
875 c sigiso(j,i) = sigiso(i,j)
878 c! IF ((wstate(k,i,j).EQ.0.0d0).AND.
879 c! & (wstate(k,j,i).NE.0.0d0)) THEN
880 c! wstate(k,i,j) = wstate(k,j,i)
881 c! ELSE IF ((wstate(k,j,i).EQ.0.0d0).AND.
882 c! & (wstate(k,i,j).NE.0.0d0)) THEN
883 c! wstate(k,j,i) = wstate(k,i,j)
888 c! THE LOOP MAKING MOMO_TABLES SYMMETRIC ENDS HERE
890 if (.not.lprint) goto 70
892 & "Parameters of the new physics-based SC-SC interaction potential"
893 write (iout,'(/7a)') 'Residues',' epsGB',' rGB',
894 & ' chi1GB',' chi2GB',' chip1GB',' chip2GB'
897 write (iout,'(2(a3,1x),1pe10.3,5(0pf10.3))')
898 & restyp(i),restyp(j),eps(i,j),sigma(i,j),chi(i,j),chi(j,i),
899 & chipp(i,j),chipp(j,i)
902 write (iout,'(/9a)') 'Residues',' alphasur1',' alphasur2',
903 & ' alphasur3',' alphasur4',' sigmap1',' sigmap2',
907 write (iout,'(2(a3,1x),8(0pf10.3))')
908 & restyp(i),restyp(j),(alphasur(k,i,j),k=1,4),
909 & sigmap(i,j),sigmap(j,i),chis(i,j),chis(j,i)
912 write (iout,'(/14a)') 'Residues',' nst ',' wst1',
913 & ' wst2',' wst3',' wst4',' dh11',' dh21',
914 & ' dh12',' dh22',' dt1',' dt2',' epsh1',
918 write (iout,'(2(a3,1x),i3,4f8.3,6f7.2,f9.5,f7.2)')
919 & restyp(i),restyp(j),nstate(i,j),(wstate(k,i,j),k=1,4),
920 & ((dhead(l,k,i,j),l=1,2),k=1,2),dtail(i,j),dtail(j,i),
921 & epshead(i,j),sig0head(i,j)
924 write (iout,'(/12a)') 'Residues',' ch1',' ch2',
925 & ' rborn1',' rborn2',' wqdip1',' wqdip2',
929 write (iout,'(2(a3,1x),2i4,5f10.3)')
930 & restyp(i),restyp(j),icharge(i),icharge(j),
931 & rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j)
934 write (iout,'(/12a)') 'Residues',
936 & ' alphpol2',' alphiso1',' alpiso2',
937 & ' alpiso3',' alpiso4',' sigiso1',' sigiso2',
941 write (iout,'(2(a3,1x),11f10.3)')
942 & restyp(i),restyp(j),alphapol(i,j),alphapol(j,i),
943 & (alphiso(k,i,j),k=1,4),sigiso(i,j),sigiso(j,i),epsintab(i,j)
948 C For the GB potential convert sigma'**2 into chi'
950 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
953 write (iout,'(/a/)') 'Parameters of the BP potential:'
954 write (iout,'(a/)') 'The epsilon array:'
955 CALL printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
956 write (iout,'(/a)') 'One-body parameters:'
957 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
959 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
960 & chip(i),alp(i),i=1,ntyp)
965 C-----------------------------------------------------------------------
966 C Calculate the "working" parameters of SC interactions.
971 alphasur(k,j,i)=alphasur(k,i,j)
972 alphiso(k,j,i)=alphiso(k,i,j)
973 wstate(k,j,i)=wstate(k,i,j)
976 wqdip(k,j,i)=wqdip(k,i,j)
980 dhead(l,k,j,i)=dhead(l,k,i,j)
988 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
989 sigma(j,i)=sigma(i,j)
990 rs0(i,j)=dwa16*sigma(i,j)
998 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
999 & 'Working parameters of the SC interactions:',
1000 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1005 IF (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4 .or. ipot.eq.6)
1015 sigeps=dsign(1.0D0,epsij)
1017 aa(i,j)=epsij*rrij*rrij
1018 bb(i,j)=-sigeps*epsij*rrij
1021 IF ((ipot.GT.2).AND.(ipot.LT.6)) THEN
1022 sigt1sq = sigma0(i)**2
1023 sigt2sq = sigma0(j)**2
1026 ratsig1 = sigt2sq/sigt1sq
1027 ratsig2 = 1.0D0/ratsig1
1028 chi(i,j) = (sigii1-1.0D0) / (sigii1+ratsig1)
1029 IF (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1030 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1034 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1035 sigmaii(i,j)=rsum_max
1036 sigmaii(j,i)=rsum_max
1038 c sigmaii(i,j)=r0(i,j)
1039 c sigmaii(j,i)=r0(i,j)
1041 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1042 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1043 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1044 augm(i,j)=epsij*r_augm**(2*expon)
1045 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1054 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1055 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1056 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1058 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3),2i3,10f8.4,
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),
1062 & icharge(i),icharge(j),chipp(i,j),chipp(j,i),
1063 & (alphasur(k,i,j),k=1,4),sigmap(i,j),sigmap(j,i),
1064 & chis(i,j),chis(j,i),
1065 & nstate(i,j),(wstate(k,i,j),k=1,4),
1066 & ((dhead(l,k,i,j),l=1,2),k=1,2),dtail(i,j),dtail(j,i),
1067 & epshead(i,j),sig0head(i,j),
1068 & rborn(i,j),(wqdip(k,i,j),k=1,2),wquad(i,j),
1069 & alphapol(i,j),alphapol(j,i),
1070 & (alphiso(k,i,j),k=1,4),sigiso(i,j)
1081 C Define the SC-p interaction constants (hard-coded; old style)
1084 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1086 c aad(i,1)=0.3D0*4.0D0**12
1087 C Following line for constants currently implemented
1088 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1089 aad(i,1)=1.5D0*4.0D0**12
1090 c aad(i,1)=0.17D0*5.6D0**12
1092 C "Soft" SC-p repulsion
1094 C Following line for constants currently implemented
1095 c aad(i,1)=0.3D0*4.0D0**6
1096 C "Hard" SC-p repulsion
1097 bad(i,1)=3.0D0*4.0D0**6
1098 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1107 C 8/9/01 Read the SC-p interaction constants from file
1110 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1113 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1114 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1115 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1116 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1120 write (iout,'(/a)') "Parameters of SC-p interactions:"
1122 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1123 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1128 C Define the constants of the disulfide bridge
1132 c Old arbitrary potential - commented out.
1137 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1138 c energy surface of diethyl disulfide.
1139 c A. Liwo and U. Kozlowska, 11/24/03
1156 c! print *, "END OF PARMREAD"
1157 c! print *, "eps = ", eps(9,9), "sigma = ", sigma(9,9)
1158 c! print *, "chi1 = ", chi(9,9), "chi2 = ", chi(9,9)
1159 c! print *, "chip1 = ", chipp(9,9), "chip2 = ", chipp(9,9)
1160 c! print *, "sig1 = ", sigmap(9,9), "sig2 = ", sigmap(9,9)
1161 c! print *, "chis1 = ", chis(9,9)," chis2 = ", chis(9,9)
1162 c! print *, "END OF PARMREAD"
1167 write (iout,'(/a)') "Disulfide bridge parameters:"
1168 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1169 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1170 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1171 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1175 111 write (iout,*) "Error reading bending energy parameters."
1177 112 write (iout,*) "Error reading rotamer energy parameters."
1179 113 write (iout,*) "Error reading torsional energy parameters."
1181 114 write (iout,*) "Error reading double torsional energy parameters."
1184 & "Error reading cumulant (multibody energy) parameters."
1186 116 write (iout,*) "Error reading electrostatic energy parameters."
1188 117 write (iout,*) "Error reading side chain interaction parameters."
1190 118 write (iout,*) "Error reading SCp interaction parameters."
1192 119 write (iout,*) "Error reading SCCOR parameters"
1195 call MPI_Finalize(Ierror)
1202 subroutine getenv_loc(var, val)
1203 character(*) var, val
1206 character(2000) line
1209 open (196,file='env',status='old',readonly,shared)
1211 c write(*,*)'looking for ',var
1212 10 read(196,*,err=11,end=11)line
1213 iread=index(line,var)
1214 c write(*,*)iread,' ',var,' ',line
1215 if (iread.eq.0) go to 10
1216 c write(*,*)'---> ',line
1222 iread=iread+ilen(var)+1
1223 read (line(iread:),*,err=12,end=12) val
1224 c write(*,*)'OK: ',var,' = ',val
1230 #elif (defined CRAY)
1231 integer lennam,lenval,ierror
1233 c getenv using a POSIX call, useful on the T3D
1234 c Sept 1996, comment out error check on advice of H. Pritchard
1237 if(lennam.le.0) stop '--error calling getenv--'
1238 call pxfgetenv(var,lennam,val,lenval,ierror)
1239 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1241 call getenv(var,val)