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)
283 C Read the parameters of the probability distribution/energy expression
284 C of the side chains.
287 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
291 dsc_inv(i)=1.0D0/dsc(i)
302 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
303 & ((blower(k,l,1),l=1,k),k=1,3)
305 read (irotam,*,end=112,err=112) bsc(j,i)
306 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
307 & ((blower(k,l,j),l=1,k),k=1,3)
314 akl=akl+blower(k,m,j)*blower(l,m,j)
325 write (iout,'(/a)') 'Parameters of side-chain local geometry'
330 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
331 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
332 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
333 & 'log h',(bsc(j,i),j=1,nlobi)
334 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
335 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
337 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
338 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
341 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
342 write (iout,'(a,f10.4,4(16x,f10.4))')
343 & 'Center ',(bsc(j,i),j=1,nlobi)
344 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
353 C Read scrot parameters for potentials determined from all-atom AM1 calculations
354 C added by Urszula Kozlowska 07/11/2007
357 read (irotam,*,end=112,err=112)
359 read (irotam,*,end=112,err=112)
362 read(irotam,*,end=112,err=112) sc_parmin(j,i)
371 C Read torsional parameters in old format
373 read (itorp,*,end=113,err=113) ntortyp,nterm_old
374 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
375 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
380 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
386 write (iout,'(/a/)') 'Torsional constants:'
389 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
390 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
396 C Read torsional parameters
398 read (itorp,*,end=113,err=113) ntortyp
399 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
400 itortyp(ntyp1)=ntortyp+1
401 c write (iout,*) 'ntortyp',ntortyp
404 read (itorp,*,end=113,err=113) nterm(i,j),nlor(i,j)
408 read (itorp,*,end=113,err=113) kk,v1(k,i,j),v2(k,i,j)
409 v0ij=v0ij+si*v1(k,i,j)
413 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
414 & vlor2(k,i,j),vlor3(k,i,j)
415 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
422 write (iout,'(/a/)') 'Torsional constants:'
425 write (iout,*) 'ityp',i,' jtyp',j
426 write (iout,*) 'Fourier constants'
428 write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
430 write (iout,*) 'Lorenz constants'
432 write (iout,'(3(1pe15.5))')
433 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
439 C 6/23/01 Read parameters for double torsionals
444 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
445 if (t1.ne.onelett(i) .or. t2.ne.onelett(j)
446 & .or. t3.ne.onelett(k)) then
447 write (iout,*) "Error in double torsional parameter file",
450 call MPI_Finalize(Ierror)
452 stop "Error in double torsional parameter file"
454 read (itordp,*,end=114,err=114) ntermd_1(i,j,k),
456 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k),l=1,
458 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k),l=1,
460 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k),l=1,
462 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k),l=1,
464 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k),
465 & v2c(m,l,i,j,k),v2s(l,m,i,j,k),v2s(m,l,i,j,k),
466 & m=1,l-1),l=1,ntermd_2(i,j,k))
472 write (iout,*) 'Constants for double torsionals'
476 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
477 & ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k)
479 write (iout,*) 'Single angles:'
480 do l=1,ntermd_1(i,j,k)
481 write (iout,'(i5,2f10.5,5x,2f10.5)') l,
482 & v1c(1,l,i,j,k),v1s(1,l,i,j,k),
483 & v1c(2,l,i,j,k),v1s(2,l,i,j,k)
486 write (iout,*) 'Pairs of angles:'
487 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
488 do l=1,ntermd_2(i,j,k)
489 write (iout,'(i5,20f10.5)')
490 & l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k))
493 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
494 do l=1,ntermd_2(i,j,k)
495 write (iout,'(i5,20f10.5)')
496 & l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k))
504 C Read of Side-chain backbone correlation parameters
505 C Modified 11 May 2012 by Adasko
508 C 5/21/07 (AL) Read coefficients of the backbone-local sidechain-local
510 read (isccor,*,end=119,err=119) nsccortyp
512 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
514 isccortyp(i)=-isccortyp(-i)
516 iscprol=isccortyp(20)
517 c write (iout,*) 'ntortyp',ntortyp
519 cc maxinter is maximum interaction sites
523 read (isccor,*,end=119,err=119)
524 &nterm_sccor(i,j),nlor_sccor(i,j)
530 nterm_sccor(-i,j)=nterm_sccor(i,j)
531 nterm_sccor(-i,-j)=nterm_sccor(i,j)
532 nterm_sccor(i,-j)=nterm_sccor(i,j)
533 do k=1,nterm_sccor(i,j)
534 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
536 if (j.eq.iscprol) then
537 if (i.eq.isccortyp(10)) then
538 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
539 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
541 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
542 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
543 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
544 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
545 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
546 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
547 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
548 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
551 if (i.eq.isccortyp(10)) then
552 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
553 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
555 if (j.eq.isccortyp(10)) then
556 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
557 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
559 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
560 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
561 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
562 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
563 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
564 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
568 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
569 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
570 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
571 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
574 do k=1,nlor_sccor(i,j)
575 read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
576 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
577 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
578 &(1+vlor3sccor(k,i,j)**2)
580 v0sccor(l,i,j)=v0ijsccor
581 v0sccor(l,-i,j)=v0ijsccor1
582 v0sccor(l,i,-j)=v0ijsccor2
583 v0sccor(l,-i,-j)=v0ijsccor3
589 read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
590 c write (iout,*) 'ntortyp',ntortyp
592 cc maxinter is maximum interaction sites
596 read (isccor,*,end=113,err=113)
597 & nterm_sccor(i,j),nlor_sccor(i,j)
601 do k=1,nterm_sccor(i,j)
602 read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j)
604 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
607 do k=1,nlor_sccor(i,j)
608 read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
609 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
610 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
611 &(1+vlor3sccor(k,i,j)**2)
613 v0sccor(l,i,j)=v0ijsccor
621 write (iout,'(/a/)') 'Torsional constants:'
624 write (iout,*) 'ityp',i,' jtyp',j
625 write (iout,*) 'Fourier constants'
626 do k=1,nterm_sccor(i,j)
627 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
629 write (iout,*) 'Lorenz constants'
630 do k=1,nlor_sccor(i,j)
631 write (iout,'(3(1pe15.5))')
632 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
639 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
640 C interaction energy of the Gly, Ala, and Pro prototypes.
644 write (iout,*) "Coefficients of the cumulants"
646 read (ifourier,*) nloctyp
648 read (ifourier,*,end=115,err=115)
649 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
651 write (iout,*) 'Type',i
652 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
678 c Ctilde(1,1,i)=0.0d0
679 c Ctilde(1,2,i)=0.0d0
680 c Ctilde(2,1,i)=0.0d0
681 c Ctilde(2,2,i)=0.0d0
694 c Dtilde(1,1,i)=0.0d0
695 c Dtilde(1,2,i)=0.0d0
696 c Dtilde(2,1,i)=0.0d0
697 c Dtilde(2,2,i)=0.0d0
698 EE(1,1,i)= b(10)+b(11)
699 EE(2,2,i)=-b(10)+b(11)
700 EE(2,1,i)= b(12)-b(13)
701 EE(1,2,i)= b(12)+b(13)
706 c ee(2,1,i)=ee(1,2,i)
708 C write(iout,*) "parm", B1(1,nloctyp+1),B1(2,nloctyp+1)
711 write (iout,*) 'Type',i
713 write(iout,*) B1(1,i),B1(2,i)
715 write(iout,*) B2(1,i),B2(2,i)
718 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
722 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
726 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
731 C Read electrostatic-interaction parameters
735 write (iout,'(/a)') 'Electrostatic interaction constants:'
736 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
737 & 'IT','JT','APP','BPP','AEL6','AEL3'
739 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
740 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
741 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
742 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
747 app (i,j)=epp(i,j)*rri*rri
748 bpp (i,j)=-2.0D0*epp(i,j)*rri
749 ael6(i,j)=elpp6(i,j)*4.2D0**6
750 ael3(i,j)=elpp3(i,j)*4.2D0**3
751 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
752 & ael6(i,j),ael3(i,j)
756 C Read side-chain interaction parameters.
758 read (isidep,*,end=117,err=117) ipot,expon
759 if (ipot.lt.1 .or. ipot.gt.5) then
760 write (iout,'(2a)') 'Error while reading SC interaction',
761 & 'potential file - unknown potential type.'
763 call MPI_Finalize(Ierror)
769 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
770 & ', exponents are ',expon,2*expon
771 goto (10,20,30,30,40) ipot
772 C----------------------- LJ potential ---------------------------------
773 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
774 & (sigma0(i),i=1,ntyp)
776 write (iout,'(/a/)') 'Parameters of the LJ potential:'
777 write (iout,'(a/)') 'The epsilon array:'
778 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
779 write (iout,'(/a)') 'One-body parameters:'
780 write (iout,'(a,4x,a)') 'residue','sigma'
781 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
784 C----------------------- LJK potential --------------------------------
785 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
786 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
788 write (iout,'(/a/)') 'Parameters of the LJK potential:'
789 write (iout,'(a/)') 'The epsilon array:'
790 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
791 write (iout,'(/a)') 'One-body parameters:'
792 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
793 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
797 C---------------------- GB or BP potential -----------------------------
798 30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
799 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
801 C For the GB potential convert sigma'**2 into chi'
804 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
808 write (iout,'(/a/)') 'Parameters of the BP potential:'
809 write (iout,'(a/)') 'The epsilon array:'
810 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
811 write (iout,'(/a)') 'One-body parameters:'
812 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
814 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
815 & chip(i),alp(i),i=1,ntyp)
818 C--------------------- GBV potential -----------------------------------
819 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
820 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
821 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
823 write (iout,'(/a/)') 'Parameters of the GBV potential:'
824 write (iout,'(a/)') 'The epsilon array:'
825 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
826 write (iout,'(/a)') 'One-body parameters:'
827 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
828 & 's||/s_|_^2',' chip ',' alph '
829 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
830 & sigii(i),chip(i),alp(i),i=1,ntyp)
834 C-----------------------------------------------------------------------
835 C Calculate the "working" parameters of SC interactions.
843 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
844 sigma(j,i)=sigma(i,j)
845 rs0(i,j)=dwa16*sigma(i,j)
849 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
850 & 'Working parameters of the SC interactions:',
851 & ' a ',' b ',' augm ',' sigma ',' r0 ',
856 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
865 sigeps=dsign(1.0D0,epsij)
867 aa(i,j)=epsij*rrij*rrij
868 bb(i,j)=-sigeps*epsij*rrij
876 ratsig1=sigt2sq/sigt1sq
877 ratsig2=1.0D0/ratsig1
878 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
879 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
880 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
884 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
885 sigmaii(i,j)=rsum_max
886 sigmaii(j,i)=rsum_max
888 c sigmaii(i,j)=r0(i,j)
889 c sigmaii(j,i)=r0(i,j)
891 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
892 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
893 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
894 augm(i,j)=epsij*r_augm**(2*expon)
895 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
902 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
903 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
904 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
910 C Define the SC-p interaction constants (hard-coded; old style)
913 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
915 c aad(i,1)=0.3D0*4.0D0**12
916 C Following line for constants currently implemented
917 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
918 aad(i,1)=1.5D0*4.0D0**12
919 c aad(i,1)=0.17D0*5.6D0**12
921 C "Soft" SC-p repulsion
923 C Following line for constants currently implemented
924 c aad(i,1)=0.3D0*4.0D0**6
925 C "Hard" SC-p repulsion
926 bad(i,1)=3.0D0*4.0D0**6
927 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
936 C 8/9/01 Read the SC-p interaction constants from file
939 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
942 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
943 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
944 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
945 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
949 write (iout,*) "Parameters of SC-p interactions:"
951 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
952 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
957 C Define the constants of the disulfide bridge
961 c Old arbitrary potential - commented out.
966 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
967 c energy surface of diethyl disulfide.
968 c A. Liwo and U. Kozlowska, 11/24/03
985 write (iout,'(/a)') "Disulfide bridge parameters:"
986 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
987 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
988 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
989 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
993 111 write (iout,*) "Error reading bending energy parameters."
995 112 write (iout,*) "Error reading rotamer energy parameters."
997 113 write (iout,*) "Error reading torsional energy parameters."
999 114 write (iout,*) "Error reading double torsional energy parameters."
1002 & "Error reading cumulant (multibody energy) parameters."
1004 116 write (iout,*) "Error reading electrostatic energy parameters."
1006 117 write (iout,*) "Error reading side chain interaction parameters."
1008 118 write (iout,*) "Error reading SCp interaction parameters."
1010 119 write (iout,*) "Error reading SCCOR parameters"
1013 call MPI_Finalize(Ierror)
1020 subroutine getenv_loc(var, val)
1021 character(*) var, val
1024 character(2000) line
1027 open (196,file='env',status='old',readonly,shared)
1029 c write(*,*)'looking for ',var
1030 10 read(196,*,err=11,end=11)line
1031 iread=index(line,var)
1032 c write(*,*)iread,' ',var,' ',line
1033 if (iread.eq.0) go to 10
1034 c write(*,*)'---> ',line
1040 iread=iread+ilen(var)+1
1041 read (line(iread:),*,err=12,end=12) val
1042 c write(*,*)'OK: ',var,' = ',val
1048 #elif (defined CRAY)
1049 integer lennam,lenval,ierror
1051 c getenv using a POSIX call, useful on the T3D
1052 c Sept 1996, comment out error check on advice of H. Pritchard
1055 if(lennam.le.0) stop '--error calling getenv--'
1056 call pxfgetenv(var,lennam,val,lenval,ierror)
1057 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1059 call getenv(var,val)