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_pdb,*,err=111,end=111)
395 & a0thet(i),(athet(j,i,1,1),j=1,2),
396 & (bthet(j,i,1,1),j=1,2)
397 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
398 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
399 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
403 athet(1,i,1,-1)=athet(1,i,1,1)
404 athet(2,i,1,-1)=athet(2,i,1,1)
405 bthet(1,i,1,-1)=-bthet(1,i,1,1)
406 bthet(2,i,1,-1)=-bthet(2,i,1,1)
407 athet(1,i,-1,1)=-athet(1,i,1,1)
408 athet(2,i,-1,1)=-athet(2,i,1,1)
409 bthet(1,i,-1,1)=bthet(1,i,1,1)
410 bthet(2,i,-1,1)=bthet(2,i,1,1)
414 athet(1,i,-1,-1)=athet(1,-i,1,1)
415 athet(2,i,-1,-1)=-athet(2,-i,1,1)
416 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
417 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
418 athet(1,i,-1,1)=athet(1,-i,1,1)
419 athet(2,i,-1,1)=-athet(2,-i,1,1)
420 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
421 bthet(2,i,-1,1)=bthet(2,-i,1,1)
422 athet(1,i,1,-1)=-athet(1,-i,1,1)
423 athet(2,i,1,-1)=athet(2,-i,1,1)
424 bthet(1,i,1,-1)=bthet(1,-i,1,1)
425 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
430 polthet(j,i)=polthet(j,-i)
433 gthet(j,i)=gthet(j,-i)
436 write (2,*) "End reading THETA_PDB"
442 C Read the parameters of the probability distribution/energy expression
443 C of the side chains.
446 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
450 dsc_inv(i)=1.0D0/dsc(i)
461 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
462 & ((blower(k,l,1),l=1,k),k=1,3)
463 censc(1,1,-i)=censc(1,1,i)
464 censc(2,1,-i)=censc(2,1,i)
465 censc(3,1,-i)=-censc(3,1,i)
468 read (irotam,*,end=112,err=112) bsc(j,i)
469 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
470 & ((blower(k,l,j),l=1,k),k=1,3)
471 censc(1,j,-i)=censc(1,j,i)
472 censc(2,j,-i)=censc(2,j,i)
473 censc(3,j,-i)=-censc(3,j,i)
474 C BSC is amplitude of Gaussian
481 akl=akl+blower(k,m,j)*blower(l,m,j)
485 if (((k.eq.3).and.(l.ne.3))
486 & .or.((l.eq.3).and.(k.ne.3))) then
487 gaussc(k,l,j,-i)=-akl
488 gaussc(l,k,j,-i)=-akl
500 write (iout,'(/a)') 'Parameters of side-chain local geometry'
505 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
506 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
507 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
508 & 'log h',(bsc(j,i),j=1,nlobi)
509 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
510 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
512 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
513 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
516 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
517 write (iout,'(a,f10.4,4(16x,f10.4))')
518 & 'Center ',(bsc(j,i),j=1,nlobi)
519 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
528 C Read scrot parameters for potentials determined from all-atom AM1 calculations
529 C added by Urszula Kozlowska 07/11/2007
532 read (irotam,*,end=112,err=112)
534 read (irotam,*,end=112,err=112)
537 read(irotam,*,end=112,err=112) sc_parmin(j,i)
542 C Read the parameters of the probability distribution/energy expression
543 C of the side chains.
545 write (2,*) "Start reading ROTAM_PDB"
548 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
552 dsc_inv(i)=1.0D0/dsc(i)
563 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
564 & ((blower(k,l,1),l=1,k),k=1,3)
566 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
567 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
568 & ((blower(k,l,j),l=1,k),k=1,3)
575 akl=akl+blower(k,m,j)*blower(l,m,j)
585 write (2,*) "Ending reading ROTAM_PDB"
591 C Read torsional parameters in old format
593 read (itorp,*,end=113,err=113) ntortyp,nterm_old
594 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
595 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
600 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
606 write (iout,'(/a/)') 'Torsional constants:'
609 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
610 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
616 C Read torsional parameters
618 read (itorp,*,end=113,err=113) ntortyp
619 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
622 itortyp(i)=-itortyp(-i)
624 c write (iout,*) 'ntortyp',ntortyp
626 do j=-ntortyp+1,ntortyp-1
627 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
629 nterm(-i,-j,iblock)=nterm(i,j,iblock)
630 nlor(-i,-j,iblock)=nlor(i,j,iblock)
633 do k=1,nterm(i,j,iblock)
634 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
636 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
637 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
638 v0ij=v0ij+si*v1(k,i,j,iblock)
640 c write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
641 c write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
642 c &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
644 do k=1,nlor(i,j,iblock)
645 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
646 & vlor2(k,i,j),vlor3(k,i,j)
647 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
650 v0(-i,-j,iblock)=v0ij
656 write (iout,'(/a/)') 'Torsional constants:'
659 write (iout,*) 'ityp',i,' jtyp',j
660 write (iout,*) 'Fourier constants'
661 do k=1,nterm(i,j,iblock)
662 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
665 write (iout,*) 'Lorenz constants'
666 do k=1,nlor(i,j,iblock)
667 write (iout,'(3(1pe15.5))')
668 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
674 C 6/23/01 Read parameters for double torsionals
678 do j=-ntortyp+1,ntortyp-1
679 do k=-ntortyp+1,ntortyp-1
680 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
681 c write (iout,*) "OK onelett",
684 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
685 & .or. t3.ne.toronelet(k)) then
686 write (iout,*) "Error in double torsional parameter file",
689 call MPI_Finalize(Ierror)
691 stop "Error in double torsional parameter file"
693 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
694 & ntermd_2(i,j,k,iblock)
695 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
696 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
697 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
698 & ntermd_1(i,j,k,iblock))
699 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
700 & ntermd_1(i,j,k,iblock))
701 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
702 & ntermd_1(i,j,k,iblock))
703 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
704 & ntermd_1(i,j,k,iblock))
705 C Martix of D parameters for one dimesional foureir series
706 do l=1,ntermd_1(i,j,k,iblock)
707 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
708 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
709 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
710 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
711 c write(iout,*) "whcodze" ,
712 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
714 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
715 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
716 & v2s(m,l,i,j,k,iblock),
717 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
718 C Martix of D parameters for two dimesional fourier series
719 do l=1,ntermd_2(i,j,k,iblock)
721 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
722 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
723 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
724 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
733 write (iout,*) 'Constants for double torsionals'
736 do j=-ntortyp+1,ntortyp-1
737 do k=-ntortyp+1,ntortyp-1
738 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
739 & ' nsingle',ntermd_1(i,j,k,iblock),
740 & ' ndouble',ntermd_2(i,j,k,iblock)
742 write (iout,*) 'Single angles:'
743 do l=1,ntermd_1(i,j,k,iblock)
744 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
745 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
746 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
747 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
750 write (iout,*) 'Pairs of angles:'
751 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
752 do l=1,ntermd_2(i,j,k,iblock)
753 write (iout,'(i5,20f10.5)')
754 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
757 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
758 do l=1,ntermd_2(i,j,k,iblock)
759 write (iout,'(i5,20f10.5)')
760 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
761 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
770 C Read of Side-chain backbone correlation parameters
771 C Modified 11 May 2012 by Adasko
774 read (isccor,*,end=119,err=119) nsccortyp
776 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
778 isccortyp(i)=-isccortyp(-i)
780 iscprol=isccortyp(20)
781 c write (iout,*) 'ntortyp',ntortyp
783 cc maxinter is maximum interaction sites
787 read (isccor,*,end=119,err=119) nterm_sccor(i,j),nlor_sccor(i,j)
793 nterm_sccor(-i,j)=nterm_sccor(i,j)
794 nterm_sccor(-i,-j)=nterm_sccor(i,j)
795 nterm_sccor(i,-j)=nterm_sccor(i,j)
796 do k=1,nterm_sccor(i,j)
797 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
799 if (j.eq.iscprol) then
800 if (i.eq.isccortyp(10)) then
801 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
802 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
804 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
805 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
806 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
807 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
808 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
809 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
810 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
811 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
814 if (i.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 if (j.eq.isccortyp(10)) then
819 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
820 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)
824 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
825 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
826 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
827 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
831 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
832 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
833 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
834 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
837 do k=1,nlor_sccor(i,j)
838 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
839 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
840 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
841 &(1+vlor3sccor(k,i,j)**2)
843 v0sccor(l,i,j)=v0ijsccor
844 v0sccor(l,-i,j)=v0ijsccor1
845 v0sccor(l,i,-j)=v0ijsccor2
846 v0sccor(l,-i,-j)=v0ijsccor3
852 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
853 c write (iout,*) 'ntortyp',ntortyp
855 cc maxinter is maximum interaction sites
859 read (isccor,*,end=119,err=119)
860 & nterm_sccor(i,j),nlor_sccor(i,j)
864 do k=1,nterm_sccor(i,j)
865 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
867 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
870 do k=1,nlor_sccor(i,j)
871 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
872 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
873 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
874 &(1+vlor3sccor(k,i,j)**2)
876 v0sccor(i,j,iblock)=v0ijsccor
884 write (iout,'(/a/)') 'Torsional constants:'
887 write (iout,*) 'ityp',i,' jtyp',j
888 write (iout,*) 'Fourier constants'
889 do k=1,nterm_sccor(i,j)
890 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
892 write (iout,*) 'Lorenz constants'
893 do k=1,nlor_sccor(i,j)
894 write (iout,'(3(1pe15.5))')
895 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
902 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
903 C interaction energy of the Gly, Ala, and Pro prototypes.
907 write (iout,*) "Coefficients of the cumulants"
909 read (ifourier,*) nloctyp
911 read (ifourier,*,end=115,err=115)
912 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
914 write (iout,*) 'Type',i
915 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
957 c Ctilde(1,1,i)=0.0d0
958 c Ctilde(1,2,i)=0.0d0
959 c Ctilde(2,1,i)=0.0d0
960 c Ctilde(2,2,i)=0.0d0
982 c Dtilde(1,1,i)=0.0d0
983 c Dtilde(1,2,i)=0.0d0
984 c Dtilde(2,1,i)=0.0d0
985 c Dtilde(2,2,i)=0.0d0
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)
990 EE(1,1,-i)= b(10)+b(11)
991 EE(2,2,-i)=-b(10)+b(11)
992 EE(2,1,-i)=-b(12)+b(13)
993 EE(1,2,-i)=-b(12)-b(13)
999 c ee(2,1,i)=ee(1,2,i)
1003 write (iout,*) 'Type',i
1005 write(iout,*) B1(1,i),B1(2,i)
1007 write(iout,*) B2(1,i),B2(2,i)
1010 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1014 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1018 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1023 C Read electrostatic-interaction parameters
1027 write (iout,'(/a)') 'Electrostatic interaction constants:'
1028 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1029 & 'IT','JT','APP','BPP','AEL6','AEL3'
1031 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1032 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1033 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1034 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1039 app (i,j)=epp(i,j)*rri*rri
1040 bpp (i,j)=-2.0D0*epp(i,j)*rri
1041 ael6(i,j)=elpp6(i,j)*4.2D0**6
1042 ael3(i,j)=elpp3(i,j)*4.2D0**3
1043 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1044 & ael6(i,j),ael3(i,j)
1048 C Read side-chain interaction parameters.
1050 read (isidep,*,end=117,err=117) ipot,expon
1051 if (ipot.lt.1 .or. ipot.gt.5) then
1052 write (iout,'(2a)') 'Error while reading SC interaction',
1053 & 'potential file - unknown potential type.'
1055 call MPI_Finalize(Ierror)
1061 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1062 & ', exponents are ',expon,2*expon
1063 goto (10,20,30,30,40) ipot
1064 C----------------------- LJ potential ---------------------------------
1065 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1066 & (sigma0(i),i=1,ntyp)
1068 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1069 write (iout,'(a/)') 'The epsilon array:'
1070 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1071 write (iout,'(/a)') 'One-body parameters:'
1072 write (iout,'(a,4x,a)') 'residue','sigma'
1073 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1076 C----------------------- LJK potential --------------------------------
1077 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1078 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1080 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1081 write (iout,'(a/)') 'The epsilon array:'
1082 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1083 write (iout,'(/a)') 'One-body parameters:'
1084 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1085 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1089 C---------------------- GB or BP potential -----------------------------
1091 read (isidep,*,end=116,err=116)(eps(i,j),j=i,ntyp)
1093 read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
1094 read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
1095 read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
1096 read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
1098 c 30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1099 c & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
1100 c & (alp(i),i=1,ntyp)
1101 C For the GB potential convert sigma'**2 into chi'
1104 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1108 write (iout,'(/a/)') 'Parameters of the BP potential:'
1109 write (iout,'(a/)') 'The epsilon array:'
1110 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1111 write (iout,'(/a)') 'One-body parameters:'
1112 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1114 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1115 & chip(i),alp(i),i=1,ntyp)
1118 C--------------------- GBV potential -----------------------------------
1119 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1120 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1121 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1123 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1124 write (iout,'(a/)') 'The epsilon array:'
1125 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1126 write (iout,'(/a)') 'One-body parameters:'
1127 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1128 & 's||/s_|_^2',' chip ',' alph '
1129 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1130 & sigii(i),chip(i),alp(i),i=1,ntyp)
1134 C-----------------------------------------------------------------------
1135 C Calculate the "working" parameters of SC interactions.
1143 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1144 sigma(j,i)=sigma(i,j)
1145 rs0(i,j)=dwa16*sigma(i,j)
1149 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1150 & 'Working parameters of the SC interactions:',
1151 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1156 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1165 sigeps=dsign(1.0D0,epsij)
1167 aa(i,j)=epsij*rrij*rrij
1168 bb(i,j)=-sigeps*epsij*rrij
1172 sigt1sq=sigma0(i)**2
1173 sigt2sq=sigma0(j)**2
1176 ratsig1=sigt2sq/sigt1sq
1177 ratsig2=1.0D0/ratsig1
1178 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1179 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1180 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1184 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1185 sigmaii(i,j)=rsum_max
1186 sigmaii(j,i)=rsum_max
1188 c sigmaii(i,j)=r0(i,j)
1189 c sigmaii(j,i)=r0(i,j)
1191 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1192 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1193 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1194 augm(i,j)=epsij*r_augm**(2*expon)
1195 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1202 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1203 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1204 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1210 C Define the SC-p interaction constants (hard-coded; old style)
1213 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1215 c aad(i,1)=0.3D0*4.0D0**12
1216 C Following line for constants currently implemented
1217 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1218 aad(i,1)=1.5D0*4.0D0**12
1219 c aad(i,1)=0.17D0*5.6D0**12
1221 C "Soft" SC-p repulsion
1223 C Following line for constants currently implemented
1224 c aad(i,1)=0.3D0*4.0D0**6
1225 C "Hard" SC-p repulsion
1226 bad(i,1)=3.0D0*4.0D0**6
1227 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1236 C 8/9/01 Read the SC-p interaction constants from file
1239 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1242 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1243 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1244 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1245 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1249 write (iout,*) "Parameters of SC-p interactions:"
1251 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1252 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1257 C Define the constants of the disulfide bridge
1261 c Old arbitrary potential - commented out.
1266 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1267 c energy surface of diethyl disulfide.
1268 c A. Liwo and U. Kozlowska, 11/24/03
1285 write (iout,'(/a)') "Disulfide bridge parameters:"
1286 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1287 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1288 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1289 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1293 111 write (iout,*) "Error reading bending energy parameters."
1295 112 write (iout,*) "Error reading rotamer energy parameters."
1297 113 write (iout,*) "Error reading torsional energy parameters."
1299 114 write (iout,*) "Error reading double torsional energy parameters."
1302 & "Error reading cumulant (multibody energy) parameters."
1304 116 write (iout,*) "Error reading electrostatic energy parameters."
1306 117 write (iout,*) "Error reading side chain interaction parameters."
1308 118 write (iout,*) "Error reading SCp interaction parameters."
1310 119 write (iout,*) "Error reading SCCOR parameters"
1313 call MPI_Finalize(Ierror)
1320 subroutine getenv_loc(var, val)
1321 character(*) var, val
1324 character(2000) line
1327 open (196,file='env',status='old',readonly,shared)
1329 c write(*,*)'looking for ',var
1330 10 read(196,*,err=11,end=11)line
1331 iread=index(line,var)
1332 c write(*,*)iread,' ',var,' ',line
1333 if (iread.eq.0) go to 10
1334 c write(*,*)'---> ',line
1340 iread=iread+ilen(var)+1
1341 read (line(iread:),*,err=12,end=12) val
1342 c write(*,*)'OK: ',var,' = ',val
1348 #elif (defined CRAY)
1349 integer lennam,lenval,ierror
1351 c getenv using a POSIX call, useful on the T3D
1352 c Sept 1996, comment out error check on advice of H. Pritchard
1355 if(lennam.le.0) stop '--error calling getenv--'
1356 call pxfgetenv(var,lennam,val,lenval,ierror)
1357 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1359 call getenv(var,val)