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"/
31 character*1 toronelet(-2:2) /"p","a","G","A","P"/
33 dimension blower(3,3,maxlob)
35 character*3 lancuch,ucase
37 C For printing parameters after they are read set the following in the UNRES
40 C setenv PRINT_PARM YES
42 C To print parameters in LaTeX format rather than as ASCII tables:
46 call getenv_loc("PRINT_PARM",lancuch)
47 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
48 call getenv_loc("LATEX",lancuch)
49 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
51 dwa16=2.0d0**(1.0d0/6.0d0)
53 C Assign virtual-bond length
58 c Read the virtual-bond parameters, masses, and moments of inertia
59 c and Stokes' radii of the peptide group and side chains
62 read (ibond,*) vbldp0,akp,mp,ip,pstok
65 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
70 dsc_inv(i)=1.0D0/dsc(i)
74 read (ibond,*) junk,vbldp0,akp,rjunk,mp,ip,pstok
76 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
77 & j=1,nbondterm(i)),msc(i),isc(i),restok(i)
82 dsc_inv(i)=1.0D0/dsc(i)
87 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
88 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',
90 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
92 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
93 & vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
95 write (iout,'(13x,3f10.5)')
96 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
102 C Read the parameters of the probability distribution/energy expression
103 C of the virtual-bond valence angles theta
106 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),
107 & (bthet(j,i,1,1),j=1,2)
108 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
109 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
110 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
114 athet(1,i,1,-1)=athet(1,i,1,1)
115 athet(2,i,1,-1)=athet(2,i,1,1)
116 bthet(1,i,1,-1)=-bthet(1,i,1,1)
117 bthet(2,i,1,-1)=-bthet(2,i,1,1)
118 athet(1,i,-1,1)=-athet(1,i,1,1)
119 athet(2,i,-1,1)=-athet(2,i,1,1)
120 bthet(1,i,-1,1)=bthet(1,i,1,1)
121 bthet(2,i,-1,1)=bthet(2,i,1,1)
125 athet(1,i,-1,-1)=athet(1,-i,1,1)
126 athet(2,i,-1,-1)=-athet(2,-i,1,1)
127 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
128 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
129 athet(1,i,-1,1)=athet(1,-i,1,1)
130 athet(2,i,-1,1)=-athet(2,-i,1,1)
131 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
132 bthet(2,i,-1,1)=bthet(2,-i,1,1)
133 athet(1,i,1,-1)=-athet(1,-i,1,1)
134 athet(2,i,1,-1)=athet(2,-i,1,1)
135 bthet(1,i,1,-1)=bthet(1,-i,1,1)
136 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
141 polthet(j,i)=polthet(j,-i)
144 gthet(j,i)=gthet(j,-i)
151 & 'Parameters of the virtual-bond valence angles:'
152 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
153 & ' ATHETA0 ',' A1 ',' A2 ',
156 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
157 & a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
159 write (iout,'(/a/9x,5a/79(1h-))')
160 & 'Parameters of the expression for sigma(theta_c):',
161 & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
162 & ' ALPH3 ',' SIGMA0C '
164 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
165 & (polthet(j,i),j=0,3),sigc0(i)
167 write (iout,'(/a/9x,5a/79(1h-))')
168 & 'Parameters of the second gaussian:',
169 & ' THETA0 ',' SIGMA0 ',' G1 ',
172 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
173 & sig0(i),(gthet(j,i),j=1,3)
177 & 'Parameters of the virtual-bond valence angles:'
178 write (iout,'(/a/9x,5a/79(1h-))')
179 & 'Coefficients of expansion',
180 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
181 & ' b1*10^1 ',' b2*10^1 '
183 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
184 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
185 & (10*bthet(j,i,1,1),j=1,2)
187 write (iout,'(/a/9x,5a/79(1h-))')
188 & 'Parameters of the expression for sigma(theta_c):',
189 & ' alpha0 ',' alph1 ',' alph2 ',
190 & ' alhp3 ',' sigma0c '
192 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
193 & (polthet(j,i),j=0,3),sigc0(i)
195 write (iout,'(/a/9x,5a/79(1h-))')
196 & 'Parameters of the second gaussian:',
197 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
200 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
201 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
207 C Read the parameters of Utheta determined from ab initio surfaces
208 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
210 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,
211 & ntheterm3,nsingle,ndouble
212 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
213 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
215 ithetyp(i)=-ithetyp(-i)
218 do i=-maxthetyp,maxthetyp
219 do j=-maxthetyp,maxthetyp
220 do k=-maxthetyp,maxthetyp
221 aa0thet(i,j,k,iblock)=0.0d0
223 aathet(l,i,j,k,iblock)=0.0d0
227 bbthet(m,l,i,j,k,iblock)=0.0d0
228 ccthet(m,l,i,j,k,iblock)=0.0d0
229 ddthet(m,l,i,j,k,iblock)=0.0d0
230 eethet(m,l,i,j,k,iblock)=0.0d0
236 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
237 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
245 c VAR:iblock means terminally blocking group 1=non-proline 2=proline
247 c VAR:ntethtyp is type of theta potentials type currently 0=glycine
248 c VAR:1=non-glicyne non-proline 2=proline
249 c VAR:negative values for D-aminoacid
251 do j=-nthetyp,nthetyp
252 do k=-nthetyp,nthetyp
253 read (ithep,'(6a)',end=111,err=111) res1
254 c VAR: aa0thet is variable describing the average value of Foureir
255 c VAR: expansion series
256 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
257 c VAR: aathet is foureir expansion in theta/2 angle for full formula
258 c VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
259 Condens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
260 read (ithep,*,end=111,err=111)
261 &(aathet(l,i,j,k,iblock),l=1,ntheterm)
262 read (ithep,*,end=111,err=111)
263 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
264 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
265 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
266 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),
268 read (ithep,*,end=111,err=111)
269 & (((ffthet(llll,lll,ll,i,j,k,iblock),
270 & ffthet(lll,llll,ll,i,j,k,iblock),
271 & ggthet(llll,lll,ll,i,j,k,iblock),
272 & ggthet(lll,llll,ll,i,j,k,iblock),
273 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
280 C For dummy ends assign glycine-type coefficients of theta-only terms; the
281 C coefficients of theta-and-gamma-dependent terms are zero.
282 C IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
283 C RECOMENTDED AFTER VERSION 3.3)
287 c aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
288 c aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
290 c aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
291 c aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
294 c aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
296 c aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
299 C AND COMMENT THE LOOPS BELOW
303 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
304 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
306 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
307 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
310 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
312 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
316 C Substitution for D aminoacids from symmetry.
319 do j=-nthetyp,nthetyp
320 do k=-nthetyp,nthetyp
321 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
323 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
327 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
328 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
329 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
330 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
336 ffthet(llll,lll,ll,i,j,k,iblock)=
337 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
338 ffthet(lll,llll,ll,i,j,k,iblock)=
339 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
340 ggthet(llll,lll,ll,i,j,k,iblock)=
341 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
342 ggthet(lll,llll,ll,i,j,k,iblock)=
343 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
352 C Control printout of the coefficients of virtual-bond-angle potentials
355 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
359 write (iout,'(//4a)')
360 & 'Type ',onelett(i),onelett(j),onelett(k)
361 write (iout,'(//a,10x,a)') " l","a[l]"
362 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
363 write (iout,'(i2,1pe15.5)')
364 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
366 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))')
367 & "b",l,"c",l,"d",l,"e",l
369 write (iout,'(i2,4(1pe15.5))') m,
370 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
371 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
375 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
376 & "f+",l,"f-",l,"g+",l,"g-",l
379 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
380 & ffthet(n,m,l,i,j,k,iblock),
381 & ffthet(m,n,l,i,j,k,iblock),
382 & ggthet(n,m,l,i,j,k,iblock),
383 & ggthet(m,n,l,i,j,k,iblock)
392 write (2,*) "Start reading THETA_PDB"
394 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),
395 & (bthet(j,i,1,1),j=1,2)
396 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
397 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
398 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
402 athet(1,i,1,-1)=athet(1,i,1,1)
403 athet(2,i,1,-1)=athet(2,i,1,1)
404 bthet(1,i,1,-1)=-bthet(1,i,1,1)
405 bthet(2,i,1,-1)=-bthet(2,i,1,1)
406 athet(1,i,-1,1)=-athet(1,i,1,1)
407 athet(2,i,-1,1)=-athet(2,i,1,1)
408 bthet(1,i,-1,1)=bthet(1,i,1,1)
409 bthet(2,i,-1,1)=bthet(2,i,1,1)
413 athet(1,i,-1,-1)=athet(1,-i,1,1)
414 athet(2,i,-1,-1)=-athet(2,-i,1,1)
415 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
416 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
417 athet(1,i,-1,1)=athet(1,-i,1,1)
418 athet(2,i,-1,1)=-athet(2,-i,1,1)
419 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
420 bthet(2,i,-1,1)=bthet(2,-i,1,1)
421 athet(1,i,1,-1)=-athet(1,-i,1,1)
422 athet(2,i,1,-1)=athet(2,-i,1,1)
423 bthet(1,i,1,-1)=bthet(1,-i,1,1)
424 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
429 polthet(j,i)=polthet(j,-i)
432 gthet(j,i)=gthet(j,-i)
435 write (2,*) "End reading THETA_PDB"
441 C Read the parameters of the probability distribution/energy expression
442 C of the side chains.
445 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
449 dsc_inv(i)=1.0D0/dsc(i)
460 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
461 & ((blower(k,l,1),l=1,k),k=1,3)
462 censc(1,1,-i)=censc(1,1,i)
463 censc(2,1,-i)=censc(2,1,i)
464 censc(3,1,-i)=-censc(3,1,i)
467 read (irotam,*,end=112,err=112) bsc(j,i)
468 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
469 & ((blower(k,l,j),l=1,k),k=1,3)
470 censc(1,j,-i)=censc(1,j,i)
471 censc(2,j,-i)=censc(2,j,i)
472 censc(3,j,-i)=-censc(3,j,i)
473 C BSC is amplitude of Gaussian
480 akl=akl+blower(k,m,j)*blower(l,m,j)
484 if (((k.eq.3).and.(l.ne.3))
485 & .or.((l.eq.3).and.(k.ne.3))) then
486 gaussc(k,l,j,-i)=-akl
487 gaussc(l,k,j,-i)=-akl
499 write (iout,'(/a)') 'Parameters of side-chain local geometry'
504 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
505 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
506 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
507 & 'log h',(bsc(j,i),j=1,nlobi)
508 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
509 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
511 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
512 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
515 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
516 write (iout,'(a,f10.4,4(16x,f10.4))')
517 & 'Center ',(bsc(j,i),j=1,nlobi)
518 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
527 C Read scrot parameters for potentials determined from all-atom AM1 calculations
528 C added by Urszula Kozlowska 07/11/2007
531 read (irotam,*,end=112,err=112)
533 read (irotam,*,end=112,err=112)
536 read(irotam,*,end=112,err=112) sc_parmin(j,i)
541 C Read the parameters of the probability distribution/energy expression
542 C of the side chains.
545 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
549 dsc_inv(i)=1.0D0/dsc(i)
560 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
561 & ((blower(k,l,1),l=1,k),k=1,3)
563 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
564 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
565 & ((blower(k,l,j),l=1,k),k=1,3)
572 akl=akl+blower(k,m,j)*blower(l,m,j)
587 C Read torsional parameters in old format
589 read (itorp,*,end=113,err=113) ntortyp,nterm_old
590 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
591 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
596 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
602 write (iout,'(/a/)') 'Torsional constants:'
605 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
606 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
612 C Read torsional parameters
614 read (itorp,*,end=113,err=113) ntortyp
615 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
618 itortyp(i)=-itortyp(-i)
620 c write (iout,*) 'ntortyp',ntortyp
622 do j=-ntortyp+1,ntortyp-1
623 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
625 nterm(-i,-j,iblock)=nterm(i,j,iblock)
626 nlor(-i,-j,iblock)=nlor(i,j,iblock)
629 do k=1,nterm(i,j,iblock)
630 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
632 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
633 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
634 v0ij=v0ij+si*v1(k,i,j,iblock)
636 c write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
637 c write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
638 c &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
640 do k=1,nlor(i,j,iblock)
641 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
642 & vlor2(k,i,j),vlor3(k,i,j)
643 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
646 v0(-i,-j,iblock)=v0ij
652 write (iout,'(/a/)') 'Torsional constants:'
655 write (iout,*) 'ityp',i,' jtyp',j
656 write (iout,*) 'Fourier constants'
657 do k=1,nterm(i,j,iblock)
658 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
661 write (iout,*) 'Lorenz constants'
662 do k=1,nlor(i,j,iblock)
663 write (iout,'(3(1pe15.5))')
664 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
670 C 6/23/01 Read parameters for double torsionals
674 do j=-ntortyp+1,ntortyp-1
675 do k=-ntortyp+1,ntortyp-1
676 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
677 c write (iout,*) "OK onelett",
680 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
681 & .or. t3.ne.toronelet(k)) then
682 write (iout,*) "Error in double torsional parameter file",
685 call MPI_Finalize(Ierror)
687 stop "Error in double torsional parameter file"
689 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
690 & ntermd_2(i,j,k,iblock)
691 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
692 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
693 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
694 & ntermd_1(i,j,k,iblock))
695 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
696 & ntermd_1(i,j,k,iblock))
697 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
698 & ntermd_1(i,j,k,iblock))
699 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
700 & ntermd_1(i,j,k,iblock))
701 C Martix of D parameters for one dimesional foureir series
702 do l=1,ntermd_1(i,j,k,iblock)
703 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
704 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
705 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
706 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
707 c write(iout,*) "whcodze" ,
708 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
710 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
711 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
712 & v2s(m,l,i,j,k,iblock),
713 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
714 C Martix of D parameters for two dimesional fourier series
715 do l=1,ntermd_2(i,j,k,iblock)
717 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
718 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
719 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
720 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
729 write (iout,*) 'Constants for double torsionals'
732 do j=-ntortyp+1,ntortyp-1
733 do k=-ntortyp+1,ntortyp-1
734 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
735 & ' nsingle',ntermd_1(i,j,k,iblock),
736 & ' ndouble',ntermd_2(i,j,k,iblock)
738 write (iout,*) 'Single angles:'
739 do l=1,ntermd_1(i,j,k,iblock)
740 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
741 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
742 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
743 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
746 write (iout,*) 'Pairs of angles:'
747 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
748 do l=1,ntermd_2(i,j,k,iblock)
749 write (iout,'(i5,20f10.5)')
750 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
753 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
754 do l=1,ntermd_2(i,j,k,iblock)
755 write (iout,'(i5,20f10.5)')
756 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
757 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
766 C Read of Side-chain backbone correlation parameters
767 C Modified 11 May 2012 by Adasko
770 read (isccor,*,end=119,err=119) nsccortyp
772 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
774 isccortyp(i)=-isccortyp(-i)
776 iscprol=isccortyp(20)
777 c write (iout,*) 'ntortyp',ntortyp
779 cc maxinter is maximum interaction sites
783 read (isccor,*,end=119,err=119) nterm_sccor(i,j),nlor_sccor(i,j)
789 nterm_sccor(-i,j)=nterm_sccor(i,j)
790 nterm_sccor(-i,-j)=nterm_sccor(i,j)
791 nterm_sccor(i,-j)=nterm_sccor(i,j)
792 do k=1,nterm_sccor(i,j)
793 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
795 if (j.eq.iscprol) then
796 if (i.eq.isccortyp(10)) then
797 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
798 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
800 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
801 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
802 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
803 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
804 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
805 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
806 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
807 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
810 if (i.eq.isccortyp(10)) then
811 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
812 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
814 if (j.eq.isccortyp(10)) then
815 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
816 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
818 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
819 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
820 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
821 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
822 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
823 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
827 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
828 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
829 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
830 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
833 do k=1,nlor_sccor(i,j)
834 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
835 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
836 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
837 &(1+vlor3sccor(k,i,j)**2)
839 v0sccor(l,i,j)=v0ijsccor
840 v0sccor(l,-i,j)=v0ijsccor1
841 v0sccor(l,i,-j)=v0ijsccor2
842 v0sccor(l,-i,-j)=v0ijsccor3
848 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
849 c write (iout,*) 'ntortyp',ntortyp
851 cc maxinter is maximum interaction sites
855 read (isccor,*,end=119,err=119)
856 & nterm_sccor(i,j),nlor_sccor(i,j)
860 do k=1,nterm_sccor(i,j)
861 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
863 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
866 do k=1,nlor_sccor(i,j)
867 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
868 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
869 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
870 &(1+vlor3sccor(k,i,j)**2)
872 v0sccor(i,j,iblock)=v0ijsccor
880 write (iout,'(/a/)') 'Torsional constants:'
883 write (iout,*) 'ityp',i,' jtyp',j
884 write (iout,*) 'Fourier constants'
885 do k=1,nterm_sccor(i,j)
886 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
888 write (iout,*) 'Lorenz constants'
889 do k=1,nlor_sccor(i,j)
890 write (iout,'(3(1pe15.5))')
891 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
898 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
899 C interaction energy of the Gly, Ala, and Pro prototypes.
903 write (iout,*) "Coefficients of the cumulants"
905 read (ifourier,*) nloctyp
907 read (ifourier,*,end=115,err=115)
908 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
910 write (iout,*) 'Type',i
911 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
953 c Ctilde(1,1,i)=0.0d0
954 c Ctilde(1,2,i)=0.0d0
955 c Ctilde(2,1,i)=0.0d0
956 c Ctilde(2,2,i)=0.0d0
978 c Dtilde(1,1,i)=0.0d0
979 c Dtilde(1,2,i)=0.0d0
980 c Dtilde(2,1,i)=0.0d0
981 c Dtilde(2,2,i)=0.0d0
982 EE(1,1,i)= b(10)+b(11)
983 EE(2,2,i)=-b(10)+b(11)
984 EE(2,1,i)= b(12)-b(13)
985 EE(1,2,i)= b(12)+b(13)
986 EE(1,1,-i)= b(10)+b(11)
987 EE(2,2,-i)=-b(10)+b(11)
988 EE(2,1,-i)=-b(12)+b(13)
989 EE(1,2,-i)=-b(12)-b(13)
995 c ee(2,1,i)=ee(1,2,i)
999 write (iout,*) 'Type',i
1001 write(iout,*) B1(1,i),B1(2,i)
1003 write(iout,*) B2(1,i),B2(2,i)
1006 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1010 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1014 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1019 C Read electrostatic-interaction parameters
1023 write (iout,'(/a)') 'Electrostatic interaction constants:'
1024 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1025 & 'IT','JT','APP','BPP','AEL6','AEL3'
1027 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1028 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1029 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1030 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1035 app (i,j)=epp(i,j)*rri*rri
1036 bpp (i,j)=-2.0D0*epp(i,j)*rri
1037 ael6(i,j)=elpp6(i,j)*4.2D0**6
1038 ael3(i,j)=elpp3(i,j)*4.2D0**3
1039 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1040 & ael6(i,j),ael3(i,j)
1044 C Read side-chain interaction parameters.
1046 read (isidep,*,end=117,err=117) ipot,expon
1047 if (ipot.lt.1 .or. ipot.gt.5) then
1048 write (iout,'(2a)') 'Error while reading SC interaction',
1049 & 'potential file - unknown potential type.'
1051 call MPI_Finalize(Ierror)
1057 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1058 & ', exponents are ',expon,2*expon
1059 goto (10,20,30,30,40) ipot
1060 C----------------------- LJ potential ---------------------------------
1061 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1062 & (sigma0(i),i=1,ntyp)
1064 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1065 write (iout,'(a/)') 'The epsilon array:'
1066 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1067 write (iout,'(/a)') 'One-body parameters:'
1068 write (iout,'(a,4x,a)') 'residue','sigma'
1069 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1072 C----------------------- LJK potential --------------------------------
1073 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1074 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1076 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1077 write (iout,'(a/)') 'The epsilon array:'
1078 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1079 write (iout,'(/a)') 'One-body parameters:'
1080 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1081 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1085 C---------------------- GB or BP potential -----------------------------
1087 read (isidep,*,end=116,err=116)(eps(i,j),j=i,ntyp)
1089 read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
1090 read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
1091 read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
1092 read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
1094 c 30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1095 c & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
1096 c & (alp(i),i=1,ntyp)
1097 C For the GB potential convert sigma'**2 into chi'
1100 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1104 write (iout,'(/a/)') 'Parameters of the BP potential:'
1105 write (iout,'(a/)') 'The epsilon array:'
1106 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1107 write (iout,'(/a)') 'One-body parameters:'
1108 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1110 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1111 & chip(i),alp(i),i=1,ntyp)
1114 C--------------------- GBV potential -----------------------------------
1115 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1116 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1117 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1119 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1120 write (iout,'(a/)') 'The epsilon array:'
1121 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1122 write (iout,'(/a)') 'One-body parameters:'
1123 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1124 & 's||/s_|_^2',' chip ',' alph '
1125 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1126 & sigii(i),chip(i),alp(i),i=1,ntyp)
1130 C-----------------------------------------------------------------------
1131 C Calculate the "working" parameters of SC interactions.
1139 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1140 sigma(j,i)=sigma(i,j)
1141 rs0(i,j)=dwa16*sigma(i,j)
1145 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1146 & 'Working parameters of the SC interactions:',
1147 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1152 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1161 sigeps=dsign(1.0D0,epsij)
1163 aa(i,j)=epsij*rrij*rrij
1164 bb(i,j)=-sigeps*epsij*rrij
1168 sigt1sq=sigma0(i)**2
1169 sigt2sq=sigma0(j)**2
1172 ratsig1=sigt2sq/sigt1sq
1173 ratsig2=1.0D0/ratsig1
1174 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1175 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1176 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1180 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1181 sigmaii(i,j)=rsum_max
1182 sigmaii(j,i)=rsum_max
1184 c sigmaii(i,j)=r0(i,j)
1185 c sigmaii(j,i)=r0(i,j)
1187 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1188 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1189 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1190 augm(i,j)=epsij*r_augm**(2*expon)
1191 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1198 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1199 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1200 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1206 C Define the SC-p interaction constants (hard-coded; old style)
1209 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1211 c aad(i,1)=0.3D0*4.0D0**12
1212 C Following line for constants currently implemented
1213 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1214 aad(i,1)=1.5D0*4.0D0**12
1215 c aad(i,1)=0.17D0*5.6D0**12
1217 C "Soft" SC-p repulsion
1219 C Following line for constants currently implemented
1220 c aad(i,1)=0.3D0*4.0D0**6
1221 C "Hard" SC-p repulsion
1222 bad(i,1)=3.0D0*4.0D0**6
1223 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1232 C 8/9/01 Read the SC-p interaction constants from file
1235 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1238 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1239 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1240 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1241 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1245 write (iout,*) "Parameters of SC-p interactions:"
1247 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1248 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1253 C Define the constants of the disulfide bridge
1257 c Old arbitrary potential - commented out.
1262 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1263 c energy surface of diethyl disulfide.
1264 c A. Liwo and U. Kozlowska, 11/24/03
1281 write (iout,'(/a)') "Disulfide bridge parameters:"
1282 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1283 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1284 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1285 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1289 111 write (iout,*) "Error reading bending energy parameters."
1291 112 write (iout,*) "Error reading rotamer energy parameters."
1293 113 write (iout,*) "Error reading torsional energy parameters."
1295 114 write (iout,*) "Error reading double torsional energy parameters."
1298 & "Error reading cumulant (multibody energy) parameters."
1300 116 write (iout,*) "Error reading electrostatic energy parameters."
1302 117 write (iout,*) "Error reading side chain interaction parameters."
1304 118 write (iout,*) "Error reading SCp interaction parameters."
1306 119 write (iout,*) "Error reading SCCOR parameters"
1309 call MPI_Finalize(Ierror)
1316 subroutine getenv_loc(var, val)
1317 character(*) var, val
1320 character(2000) line
1323 open (196,file='env',status='old',readonly,shared)
1325 c write(*,*)'looking for ',var
1326 10 read(196,*,err=11,end=11)line
1327 iread=index(line,var)
1328 c write(*,*)iread,' ',var,' ',line
1329 if (iread.eq.0) go to 10
1330 c write(*,*)'---> ',line
1336 iread=iread+ilen(var)+1
1337 read (line(iread:),*,err=12,end=12) val
1338 c write(*,*)'OK: ',var,' = ',val
1344 #elif (defined CRAY)
1345 integer lennam,lenval,ierror
1347 c getenv using a POSIX call, useful on the T3D
1348 c Sept 1996, comment out error check on advice of H. Pritchard
1351 if(lennam.le.0) stop '--error calling getenv--'
1352 call pxfgetenv(var,lennam,val,lenval,ierror)
1353 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1355 call getenv(var,val)