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)
710 write (iout,*) 'Type',i
712 write(iout,*) B1(1,i),B1(2,i)
714 write(iout,*) B2(1,i),B2(2,i)
717 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
721 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
725 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
730 C Read electrostatic-interaction parameters
734 write (iout,'(/a)') 'Electrostatic interaction constants:'
735 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
736 & 'IT','JT','APP','BPP','AEL6','AEL3'
738 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
739 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
740 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
741 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
746 app (i,j)=epp(i,j)*rri*rri
747 bpp (i,j)=-2.0D0*epp(i,j)*rri
748 ael6(i,j)=elpp6(i,j)*4.2D0**6
749 ael3(i,j)=elpp3(i,j)*4.2D0**3
750 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
751 & ael6(i,j),ael3(i,j)
755 C Read side-chain interaction parameters.
757 read (isidep,*,end=117,err=117) ipot,expon
758 if (ipot.lt.1 .or. ipot.gt.5) then
759 write (iout,'(2a)') 'Error while reading SC interaction',
760 & 'potential file - unknown potential type.'
762 call MPI_Finalize(Ierror)
768 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
769 & ', exponents are ',expon,2*expon
770 goto (10,20,30,30,40) ipot
771 C----------------------- LJ potential ---------------------------------
772 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
773 & (sigma0(i),i=1,ntyp)
775 write (iout,'(/a/)') 'Parameters of the LJ potential:'
776 write (iout,'(a/)') 'The epsilon array:'
777 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
778 write (iout,'(/a)') 'One-body parameters:'
779 write (iout,'(a,4x,a)') 'residue','sigma'
780 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
783 C----------------------- LJK potential --------------------------------
784 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
785 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
787 write (iout,'(/a/)') 'Parameters of the LJK potential:'
788 write (iout,'(a/)') 'The epsilon array:'
789 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
790 write (iout,'(/a)') 'One-body parameters:'
791 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
792 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
796 C---------------------- GB or BP potential -----------------------------
797 30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
798 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
800 C For the GB potential convert sigma'**2 into chi'
803 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
807 write (iout,'(/a/)') 'Parameters of the BP potential:'
808 write (iout,'(a/)') 'The epsilon array:'
809 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
810 write (iout,'(/a)') 'One-body parameters:'
811 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
813 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
814 & chip(i),alp(i),i=1,ntyp)
817 C--------------------- GBV potential -----------------------------------
818 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
819 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
820 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
822 write (iout,'(/a/)') 'Parameters of the GBV potential:'
823 write (iout,'(a/)') 'The epsilon array:'
824 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
825 write (iout,'(/a)') 'One-body parameters:'
826 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
827 & 's||/s_|_^2',' chip ',' alph '
828 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
829 & sigii(i),chip(i),alp(i),i=1,ntyp)
833 C-----------------------------------------------------------------------
834 C Calculate the "working" parameters of SC interactions.
842 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
843 sigma(j,i)=sigma(i,j)
844 rs0(i,j)=dwa16*sigma(i,j)
848 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
849 & 'Working parameters of the SC interactions:',
850 & ' a ',' b ',' augm ',' sigma ',' r0 ',
855 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
864 sigeps=dsign(1.0D0,epsij)
866 aa(i,j)=epsij*rrij*rrij
867 bb(i,j)=-sigeps*epsij*rrij
875 ratsig1=sigt2sq/sigt1sq
876 ratsig2=1.0D0/ratsig1
877 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
878 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
879 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
883 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
884 sigmaii(i,j)=rsum_max
885 sigmaii(j,i)=rsum_max
887 c sigmaii(i,j)=r0(i,j)
888 c sigmaii(j,i)=r0(i,j)
890 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
891 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
892 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
893 augm(i,j)=epsij*r_augm**(2*expon)
894 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
901 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
902 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
903 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
909 C Define the SC-p interaction constants (hard-coded; old style)
912 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
914 c aad(i,1)=0.3D0*4.0D0**12
915 C Following line for constants currently implemented
916 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
917 aad(i,1)=1.5D0*4.0D0**12
918 c aad(i,1)=0.17D0*5.6D0**12
920 C "Soft" SC-p repulsion
922 C Following line for constants currently implemented
923 c aad(i,1)=0.3D0*4.0D0**6
924 C "Hard" SC-p repulsion
925 bad(i,1)=3.0D0*4.0D0**6
926 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
935 C 8/9/01 Read the SC-p interaction constants from file
938 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
941 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
942 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
943 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
944 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
948 write (iout,*) "Parameters of SC-p interactions:"
950 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
951 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
956 C Define the constants of the disulfide bridge
960 c Old arbitrary potential - commented out.
965 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
966 c energy surface of diethyl disulfide.
967 c A. Liwo and U. Kozlowska, 11/24/03
984 write (iout,'(/a)') "Disulfide bridge parameters:"
985 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
986 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
987 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
988 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
992 111 write (iout,*) "Error reading bending energy parameters."
994 112 write (iout,*) "Error reading rotamer energy parameters."
996 113 write (iout,*) "Error reading torsional energy parameters."
998 114 write (iout,*) "Error reading double torsional energy parameters."
1001 & "Error reading cumulant (multibody energy) parameters."
1003 116 write (iout,*) "Error reading electrostatic energy parameters."
1005 117 write (iout,*) "Error reading side chain interaction parameters."
1007 118 write (iout,*) "Error reading SCp interaction parameters."
1009 119 write (iout,*) "Error reading SCCOR parameters"
1012 call MPI_Finalize(Ierror)
1019 subroutine getenv_loc(var, val)
1020 character(*) var, val
1023 character(2000) line
1026 open (196,file='env',status='old',readonly,shared)
1028 c write(*,*)'looking for ',var
1029 10 read(196,*,err=11,end=11)line
1030 iread=index(line,var)
1031 c write(*,*)iread,' ',var,' ',line
1032 if (iread.eq.0) go to 10
1033 c write(*,*)'---> ',line
1039 iread=iread+ilen(var)+1
1040 read (line(iread:),*,err=12,end=12) val
1041 c write(*,*)'OK: ',var,' = ',val
1047 #elif (defined CRAY)
1048 integer lennam,lenval,ierror
1050 c getenv using a POSIX call, useful on the T3D
1051 c Sept 1996, comment out error check on advice of H. Pritchard
1054 if(lennam.le.0) stop '--error calling getenv--'
1055 call pxfgetenv(var,lennam,val,lenval,ierror)
1056 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1058 call getenv(var,val)