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 c write (iout,*) "OK onelett",
500 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
501 & .or. t3.ne.toronelet(k)) then
502 write (iout,*) "Error in double torsional parameter file",
505 call MPI_Finalize(Ierror)
507 stop "Error in double torsional parameter file"
509 read (itordp,*,end=114,err=114) ntermd_1(i,j,k),
511 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k),l=1,
513 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k),l=1,
515 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k),l=1,
517 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k),l=1,
519 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k),
520 & v2c(m,l,i,j,k),v2s(l,m,i,j,k),v2s(m,l,i,j,k),
521 & m=1,l-1),l=1,ntermd_2(i,j,k))
527 write (iout,*) 'Constants for double torsionals'
531 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
532 & ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k)
534 write (iout,*) 'Single angles:'
535 do l=1,ntermd_1(i,j,k)
536 write (iout,'(i5,2f10.5,5x,2f10.5)') l,
537 & v1c(1,l,i,j,k),v1s(1,l,i,j,k),
538 & v1c(2,l,i,j,k),v1s(2,l,i,j,k)
541 write (iout,*) 'Pairs of angles:'
542 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
543 do l=1,ntermd_2(i,j,k)
544 write (iout,'(i5,20f10.5)')
545 & l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k))
548 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
549 do l=1,ntermd_2(i,j,k)
550 write (iout,'(i5,20f10.5)')
551 & l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k))
559 C Read of Side-chain backbone correlation parameters
560 C Modified 11 May 2012 by Adasko
563 read (isccor,*,end=113,err=113) nsccortyp
564 read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
566 isccortyp(i)=-isccortyp(-i)
568 iscprol=isccortyp(20)
569 c write (iout,*) 'ntortyp',ntortyp
571 cc maxinter is maximum interaction sites
575 read (isccor,*,end=113,err=113) nterm_sccor(i,j),nlor_sccor(i,j)
581 nterm_sccor(-i,j)=nterm_sccor(i,j)
582 nterm_sccor(-i,-j)=nterm_sccor(i,j)
583 nterm_sccor(i,-j)=nterm_sccor(i,j)
584 do k=1,nterm_sccor(i,j)
585 read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j)
587 if (j.eq.iscprol) then
588 if (i.eq.isccortyp(10)) then
589 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
590 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
592 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
593 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
594 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
595 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
596 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
597 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
598 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
599 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
602 if (i.eq.isccortyp(10)) then
603 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
604 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
606 if (j.eq.isccortyp(10)) then
607 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
608 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
610 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
611 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
612 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
613 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
614 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
615 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
619 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
620 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
621 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
622 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
625 do k=1,nlor_sccor(i,j)
626 read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
627 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
628 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
629 &(1+vlor3sccor(k,i,j)**2)
631 v0sccor(l,i,j)=v0ijsccor
632 v0sccor(l,-i,j)=v0ijsccor1
633 v0sccor(l,i,-j)=v0ijsccor2
634 v0sccor(l,-i,-j)=v0ijsccor3
640 read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
641 c write (iout,*) 'ntortyp',ntortyp
643 cc maxinter is maximum interaction sites
647 read (isccor,*,end=113,err=113)
648 & nterm_sccor(i,j),nlor_sccor(i,j)
652 do k=1,nterm_sccor(i,j)
653 read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j)
655 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
658 do k=1,nlor_sccor(i,j)
659 read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
660 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
661 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
662 &(1+vlor3sccor(k,i,j)**2)
664 v0sccor(i,j)=v0ijsccor
672 write (iout,'(/a/)') 'Torsional constants:'
675 write (iout,*) 'ityp',i,' jtyp',j
676 write (iout,*) 'Fourier constants'
677 do k=1,nterm_sccor(i,j)
678 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
680 write (iout,*) 'Lorenz constants'
681 do k=1,nlor_sccor(i,j)
682 write (iout,'(3(1pe15.5))')
683 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
690 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
691 C interaction energy of the Gly, Ala, and Pro prototypes.
695 write (iout,*) "Coefficients of the cumulants"
697 read (ifourier,*) nloctyp
699 read (ifourier,*,end=115,err=115)
700 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
702 write (iout,*) 'Type',i
703 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
731 c Ctilde(1,1,i)=0.0d0
732 c Ctilde(1,2,i)=0.0d0
733 c Ctilde(2,1,i)=0.0d0
734 c Ctilde(2,2,i)=0.0d0
747 c Dtilde(1,1,i)=0.0d0
748 c Dtilde(1,2,i)=0.0d0
749 c Dtilde(2,1,i)=0.0d0
750 c Dtilde(2,2,i)=0.0d0
751 EE(1,1,i)= b(10)+b(11)
752 EE(2,2,i)=-b(10)+b(11)
753 EE(2,1,i)= b(12)-b(13)
754 EE(1,2,i)= b(12)+b(13)
759 c ee(2,1,i)=ee(1,2,i)
763 write (iout,*) 'Type',i
765 write(iout,*) B1(1,i),B1(2,i)
767 write(iout,*) B2(1,i),B2(2,i)
770 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
774 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
778 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
783 C Read electrostatic-interaction parameters
787 write (iout,'(/a)') 'Electrostatic interaction constants:'
788 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
789 & 'IT','JT','APP','BPP','AEL6','AEL3'
791 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
792 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
793 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
794 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
799 app (i,j)=epp(i,j)*rri*rri
800 bpp (i,j)=-2.0D0*epp(i,j)*rri
801 ael6(i,j)=elpp6(i,j)*4.2D0**6
802 ael3(i,j)=elpp3(i,j)*4.2D0**3
803 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
804 & ael6(i,j),ael3(i,j)
808 C Read side-chain interaction parameters.
810 read (isidep,*,end=117,err=117) ipot,expon
811 if (ipot.lt.1 .or. ipot.gt.5) then
812 write (iout,'(2a)') 'Error while reading SC interaction',
813 & 'potential file - unknown potential type.'
815 call MPI_Finalize(Ierror)
821 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
822 & ', exponents are ',expon,2*expon
823 goto (10,20,30,30,40) ipot
824 C----------------------- LJ potential ---------------------------------
825 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
826 & (sigma0(i),i=1,ntyp)
828 write (iout,'(/a/)') 'Parameters of the LJ potential:'
829 write (iout,'(a/)') 'The epsilon array:'
830 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
831 write (iout,'(/a)') 'One-body parameters:'
832 write (iout,'(a,4x,a)') 'residue','sigma'
833 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
836 C----------------------- LJK potential --------------------------------
837 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
838 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
840 write (iout,'(/a/)') 'Parameters of the LJK potential:'
841 write (iout,'(a/)') 'The epsilon array:'
842 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
843 write (iout,'(/a)') 'One-body parameters:'
844 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
845 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
849 C---------------------- GB or BP potential -----------------------------
850 30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
851 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
853 C For the GB potential convert sigma'**2 into chi'
856 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
860 write (iout,'(/a/)') 'Parameters of the BP potential:'
861 write (iout,'(a/)') 'The epsilon array:'
862 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
863 write (iout,'(/a)') 'One-body parameters:'
864 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
866 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
867 & chip(i),alp(i),i=1,ntyp)
870 C--------------------- GBV potential -----------------------------------
871 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
872 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
873 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
875 write (iout,'(/a/)') 'Parameters of the GBV potential:'
876 write (iout,'(a/)') 'The epsilon array:'
877 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
878 write (iout,'(/a)') 'One-body parameters:'
879 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
880 & 's||/s_|_^2',' chip ',' alph '
881 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
882 & sigii(i),chip(i),alp(i),i=1,ntyp)
886 C-----------------------------------------------------------------------
887 C Calculate the "working" parameters of SC interactions.
895 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
896 sigma(j,i)=sigma(i,j)
897 rs0(i,j)=dwa16*sigma(i,j)
901 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
902 & 'Working parameters of the SC interactions:',
903 & ' a ',' b ',' augm ',' sigma ',' r0 ',
908 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
917 sigeps=dsign(1.0D0,epsij)
919 aa(i,j)=epsij*rrij*rrij
920 bb(i,j)=-sigeps*epsij*rrij
928 ratsig1=sigt2sq/sigt1sq
929 ratsig2=1.0D0/ratsig1
930 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
931 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
932 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
936 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
937 sigmaii(i,j)=rsum_max
938 sigmaii(j,i)=rsum_max
940 c sigmaii(i,j)=r0(i,j)
941 c sigmaii(j,i)=r0(i,j)
943 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
944 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
945 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
946 augm(i,j)=epsij*r_augm**(2*expon)
947 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
954 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
955 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
956 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
962 C Define the SC-p interaction constants (hard-coded; old style)
965 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
967 c aad(i,1)=0.3D0*4.0D0**12
968 C Following line for constants currently implemented
969 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
970 aad(i,1)=1.5D0*4.0D0**12
971 c aad(i,1)=0.17D0*5.6D0**12
973 C "Soft" SC-p repulsion
975 C Following line for constants currently implemented
976 c aad(i,1)=0.3D0*4.0D0**6
977 C "Hard" SC-p repulsion
978 bad(i,1)=3.0D0*4.0D0**6
979 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
988 C 8/9/01 Read the SC-p interaction constants from file
991 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
994 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
995 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
996 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
997 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1001 write (iout,*) "Parameters of SC-p interactions:"
1003 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1004 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1009 C Define the constants of the disulfide bridge
1013 c Old arbitrary potential - commented out.
1018 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1019 c energy surface of diethyl disulfide.
1020 c A. Liwo and U. Kozlowska, 11/24/03
1037 write (iout,'(/a)') "Disulfide bridge parameters:"
1038 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1039 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1040 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1041 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1045 111 write (iout,*) "Error reading bending energy parameters."
1047 112 write (iout,*) "Error reading rotamer energy parameters."
1049 113 write (iout,*) "Error reading torsional energy parameters."
1051 114 write (iout,*) "Error reading double torsional energy parameters."
1054 & "Error reading cumulant (multibody energy) parameters."
1056 116 write (iout,*) "Error reading electrostatic energy parameters."
1058 117 write (iout,*) "Error reading side chain interaction parameters."
1060 118 write (iout,*) "Error reading SCp interaction parameters."
1062 119 write (iout,*) "Error reading SCCOR parameters"
1065 call MPI_Finalize(Ierror)
1072 subroutine getenv_loc(var, val)
1073 character(*) var, val
1076 character(2000) line
1079 open (196,file='env',status='old',readonly,shared)
1081 c write(*,*)'looking for ',var
1082 10 read(196,*,err=11,end=11)line
1083 iread=index(line,var)
1084 c write(*,*)iread,' ',var,' ',line
1085 if (iread.eq.0) go to 10
1086 c write(*,*)'---> ',line
1092 iread=iread+ilen(var)+1
1093 read (line(iread:),*,err=12,end=12) val
1094 c write(*,*)'OK: ',var,' = ',val
1100 #elif (defined CRAY)
1101 integer lennam,lenval,ierror
1103 c getenv using a POSIX call, useful on the T3D
1104 c Sept 1996, comment out error check on advice of H. Pritchard
1107 if(lennam.le.0) stop '--error calling getenv--'
1108 call pxfgetenv(var,lennam,val,lenval,ierror)
1109 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1111 call getenv(var,val)