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
565 read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
567 isccortyp(i)=-isccortyp(-i)
569 iscprol=isccortyp(20)
570 c write (iout,*) 'ntortyp',ntortyp
572 cc maxinter is maximum interaction sites
576 read (isccor,*,end=113,err=113) nterm_sccor(i,j),nlor_sccor(i,j)
582 nterm_sccor(-i,j)=nterm_sccor(i,j)
583 nterm_sccor(-i,-j)=nterm_sccor(i,j)
584 nterm_sccor(i,-j)=nterm_sccor(i,j)
585 do k=1,nterm_sccor(i,j)
586 read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j)
588 if (j.eq.iscprol) then
589 if (i.eq.isccortyp(10)) then
590 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
591 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
593 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
594 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
595 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
596 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
597 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
598 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
599 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
600 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
603 if (i.eq.isccortyp(10)) then
604 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
605 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
607 if (j.eq.isccortyp(10)) then
608 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
609 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
611 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
612 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
613 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
614 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
615 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
616 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
620 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
621 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
622 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
623 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
626 do k=1,nlor_sccor(i,j)
627 read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
628 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
629 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
630 &(1+vlor3sccor(k,i,j)**2)
632 v0sccor(l,i,j)=v0ijsccor
633 v0sccor(l,-i,j)=v0ijsccor1
634 v0sccor(l,i,-j)=v0ijsccor2
635 v0sccor(l,-i,-j)=v0ijsccor3
641 read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
642 c write (iout,*) 'ntortyp',ntortyp
644 cc maxinter is maximum interaction sites
648 read (isccor,*,end=113,err=113)
649 & nterm_sccor(i,j),nlor_sccor(i,j)
653 do k=1,nterm_sccor(i,j)
654 read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j)
656 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
659 do k=1,nlor_sccor(i,j)
660 read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
661 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
662 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
663 &(1+vlor3sccor(k,i,j)**2)
665 v0sccor(i,j)=v0ijsccor
673 write (iout,'(/a/)') 'Torsional constants:'
676 write (iout,*) 'ityp',i,' jtyp',j
677 write (iout,*) 'Fourier constants'
678 do k=1,nterm_sccor(i,j)
679 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
681 write (iout,*) 'Lorenz constants'
682 do k=1,nlor_sccor(i,j)
683 write (iout,'(3(1pe15.5))')
684 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
691 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
692 C interaction energy of the Gly, Ala, and Pro prototypes.
696 write (iout,*) "Coefficients of the cumulants"
698 read (ifourier,*) nloctyp
700 read (ifourier,*,end=115,err=115)
701 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
703 write (iout,*) 'Type',i
704 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
732 c Ctilde(1,1,i)=0.0d0
733 c Ctilde(1,2,i)=0.0d0
734 c Ctilde(2,1,i)=0.0d0
735 c Ctilde(2,2,i)=0.0d0
748 c Dtilde(1,1,i)=0.0d0
749 c Dtilde(1,2,i)=0.0d0
750 c Dtilde(2,1,i)=0.0d0
751 c Dtilde(2,2,i)=0.0d0
752 EE(1,1,i)= b(10)+b(11)
753 EE(2,2,i)=-b(10)+b(11)
754 EE(2,1,i)= b(12)-b(13)
755 EE(1,2,i)= b(12)+b(13)
760 c ee(2,1,i)=ee(1,2,i)
764 write (iout,*) 'Type',i
766 write(iout,*) B1(1,i),B1(2,i)
768 write(iout,*) B2(1,i),B2(2,i)
771 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
775 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
779 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
784 C Read electrostatic-interaction parameters
788 write (iout,'(/a)') 'Electrostatic interaction constants:'
789 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
790 & 'IT','JT','APP','BPP','AEL6','AEL3'
792 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
793 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
794 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
795 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
800 app (i,j)=epp(i,j)*rri*rri
801 bpp (i,j)=-2.0D0*epp(i,j)*rri
802 ael6(i,j)=elpp6(i,j)*4.2D0**6
803 ael3(i,j)=elpp3(i,j)*4.2D0**3
804 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
805 & ael6(i,j),ael3(i,j)
809 C Read side-chain interaction parameters.
811 read (isidep,*,end=117,err=117) ipot,expon
812 if (ipot.lt.1 .or. ipot.gt.5) then
813 write (iout,'(2a)') 'Error while reading SC interaction',
814 & 'potential file - unknown potential type.'
816 call MPI_Finalize(Ierror)
822 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
823 & ', exponents are ',expon,2*expon
824 goto (10,20,30,30,40) ipot
825 C----------------------- LJ potential ---------------------------------
826 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
827 & (sigma0(i),i=1,ntyp)
829 write (iout,'(/a/)') 'Parameters of the LJ potential:'
830 write (iout,'(a/)') 'The epsilon array:'
831 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
832 write (iout,'(/a)') 'One-body parameters:'
833 write (iout,'(a,4x,a)') 'residue','sigma'
834 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
837 C----------------------- LJK potential --------------------------------
838 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
839 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
841 write (iout,'(/a/)') 'Parameters of the LJK potential:'
842 write (iout,'(a/)') 'The epsilon array:'
843 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
844 write (iout,'(/a)') 'One-body parameters:'
845 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
846 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
850 C---------------------- GB or BP potential -----------------------------
851 30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
852 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
854 C For the GB potential convert sigma'**2 into chi'
857 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
861 write (iout,'(/a/)') 'Parameters of the BP potential:'
862 write (iout,'(a/)') 'The epsilon array:'
863 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
864 write (iout,'(/a)') 'One-body parameters:'
865 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
867 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
868 & chip(i),alp(i),i=1,ntyp)
871 C--------------------- GBV potential -----------------------------------
872 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
873 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
874 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
876 write (iout,'(/a/)') 'Parameters of the GBV potential:'
877 write (iout,'(a/)') 'The epsilon array:'
878 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
879 write (iout,'(/a)') 'One-body parameters:'
880 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
881 & 's||/s_|_^2',' chip ',' alph '
882 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
883 & sigii(i),chip(i),alp(i),i=1,ntyp)
887 C-----------------------------------------------------------------------
888 C Calculate the "working" parameters of SC interactions.
896 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
897 sigma(j,i)=sigma(i,j)
898 rs0(i,j)=dwa16*sigma(i,j)
902 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
903 & 'Working parameters of the SC interactions:',
904 & ' a ',' b ',' augm ',' sigma ',' r0 ',
909 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
918 sigeps=dsign(1.0D0,epsij)
920 aa(i,j)=epsij*rrij*rrij
921 bb(i,j)=-sigeps*epsij*rrij
929 ratsig1=sigt2sq/sigt1sq
930 ratsig2=1.0D0/ratsig1
931 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
932 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
933 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
937 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
938 sigmaii(i,j)=rsum_max
939 sigmaii(j,i)=rsum_max
941 c sigmaii(i,j)=r0(i,j)
942 c sigmaii(j,i)=r0(i,j)
944 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
945 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
946 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
947 augm(i,j)=epsij*r_augm**(2*expon)
948 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
955 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
956 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
957 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
963 C Define the SC-p interaction constants (hard-coded; old style)
966 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
968 c aad(i,1)=0.3D0*4.0D0**12
969 C Following line for constants currently implemented
970 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
971 aad(i,1)=1.5D0*4.0D0**12
972 c aad(i,1)=0.17D0*5.6D0**12
974 C "Soft" SC-p repulsion
976 C Following line for constants currently implemented
977 c aad(i,1)=0.3D0*4.0D0**6
978 C "Hard" SC-p repulsion
979 bad(i,1)=3.0D0*4.0D0**6
980 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
989 C 8/9/01 Read the SC-p interaction constants from file
992 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
995 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
996 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
997 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
998 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1002 write (iout,*) "Parameters of SC-p interactions:"
1004 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1005 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1010 C Define the constants of the disulfide bridge
1014 c Old arbitrary potential - commented out.
1019 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1020 c energy surface of diethyl disulfide.
1021 c A. Liwo and U. Kozlowska, 11/24/03
1038 write (iout,'(/a)') "Disulfide bridge parameters:"
1039 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1040 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1041 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1042 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1046 111 write (iout,*) "Error reading bending energy parameters."
1048 112 write (iout,*) "Error reading rotamer energy parameters."
1050 113 write (iout,*) "Error reading torsional energy parameters."
1052 114 write (iout,*) "Error reading double torsional energy parameters."
1055 & "Error reading cumulant (multibody energy) parameters."
1057 116 write (iout,*) "Error reading electrostatic energy parameters."
1059 117 write (iout,*) "Error reading side chain interaction parameters."
1061 118 write (iout,*) "Error reading SCp interaction parameters."
1063 119 write (iout,*) "Error reading SCCOR parameters"
1066 call MPI_Finalize(Ierror)
1073 subroutine getenv_loc(var, val)
1074 character(*) var, val
1077 character(2000) line
1080 open (196,file='env',status='old',readonly,shared)
1082 c write(*,*)'looking for ',var
1083 10 read(196,*,err=11,end=11)line
1084 iread=index(line,var)
1085 c write(*,*)iread,' ',var,' ',line
1086 if (iread.eq.0) go to 10
1087 c write(*,*)'---> ',line
1093 iread=iread+ilen(var)+1
1094 read (line(iread:),*,err=12,end=12) val
1095 c write(*,*)'OK: ',var,' = ',val
1101 #elif (defined CRAY)
1102 integer lennam,lenval,ierror
1104 c getenv using a POSIX call, useful on the T3D
1105 c Sept 1996, comment out error check on advice of H. Pritchard
1108 if(lennam.le.0) stop '--error calling getenv--'
1109 call pxfgetenv(var,lennam,val,lenval,ierror)
1110 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1112 call getenv(var,val)