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,vbldpdum,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,vbldpdum,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)
100 C reading lipid parameters
101 read(iliptranpar,*) pepliptran
103 read(iliptranpar,*) liptranene(i)
108 C Read the parameters of the probability distribution/energy expression
109 C of the virtual-bond valence angles theta
112 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),
113 & (bthet(j,i,1,1),j=1,2)
114 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
115 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
116 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
120 athet(1,i,1,-1)=athet(1,i,1,1)
121 athet(2,i,1,-1)=athet(2,i,1,1)
122 bthet(1,i,1,-1)=-bthet(1,i,1,1)
123 bthet(2,i,1,-1)=-bthet(2,i,1,1)
124 athet(1,i,-1,1)=-athet(1,i,1,1)
125 athet(2,i,-1,1)=-athet(2,i,1,1)
126 bthet(1,i,-1,1)=bthet(1,i,1,1)
127 bthet(2,i,-1,1)=bthet(2,i,1,1)
131 athet(1,i,-1,-1)=athet(1,-i,1,1)
132 athet(2,i,-1,-1)=-athet(2,-i,1,1)
133 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
134 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
135 athet(1,i,-1,1)=athet(1,-i,1,1)
136 athet(2,i,-1,1)=-athet(2,-i,1,1)
137 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
138 bthet(2,i,-1,1)=bthet(2,-i,1,1)
139 athet(1,i,1,-1)=-athet(1,-i,1,1)
140 athet(2,i,1,-1)=athet(2,-i,1,1)
141 bthet(1,i,1,-1)=bthet(1,-i,1,1)
142 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
147 polthet(j,i)=polthet(j,-i)
150 gthet(j,i)=gthet(j,-i)
158 & 'Parameters of the virtual-bond valence angles:'
159 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
160 & ' ATHETA0 ',' A1 ',' A2 ',
163 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
164 & a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
166 write (iout,'(/a/9x,5a/79(1h-))')
167 & 'Parameters of the expression for sigma(theta_c):',
168 & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
169 & ' ALPH3 ',' SIGMA0C '
171 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
172 & (polthet(j,i),j=0,3),sigc0(i)
174 write (iout,'(/a/9x,5a/79(1h-))')
175 & 'Parameters of the second gaussian:',
176 & ' THETA0 ',' SIGMA0 ',' G1 ',
179 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
180 & sig0(i),(gthet(j,i),j=1,3)
184 & 'Parameters of the virtual-bond valence angles:'
185 write (iout,'(/a/9x,5a/79(1h-))')
186 & 'Coefficients of expansion',
187 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
188 & ' b1*10^1 ',' b2*10^1 '
190 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
191 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
192 & (10*bthet(j,i,1,1),j=1,2)
194 write (iout,'(/a/9x,5a/79(1h-))')
195 & 'Parameters of the expression for sigma(theta_c):',
196 & ' alpha0 ',' alph1 ',' alph2 ',
197 & ' alhp3 ',' sigma0c '
199 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
200 & (polthet(j,i),j=0,3),sigc0(i)
202 write (iout,'(/a/9x,5a/79(1h-))')
203 & 'Parameters of the second gaussian:',
204 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
207 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
208 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
214 C Read the parameters of Utheta determined from ab initio surfaces
215 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
217 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,
218 & ntheterm3,nsingle,ndouble
219 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
220 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
222 ithetyp(i)=-ithetyp(-i)
225 do i=-maxthetyp,maxthetyp
226 do j=-maxthetyp,maxthetyp
227 do k=-maxthetyp,maxthetyp
228 aa0thet(i,j,k,iblock)=0.0d0
230 aathet(l,i,j,k,iblock)=0.0d0
234 bbthet(m,l,i,j,k,iblock)=0.0d0
235 ccthet(m,l,i,j,k,iblock)=0.0d0
236 ddthet(m,l,i,j,k,iblock)=0.0d0
237 eethet(m,l,i,j,k,iblock)=0.0d0
243 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
244 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
252 c VAR:iblock means terminally blocking group 1=non-proline 2=proline
254 c VAR:ntethtyp is type of theta potentials type currently 0=glycine
255 c VAR:1=non-glicyne non-proline 2=proline
256 c VAR:negative values for D-aminoacid
258 do j=-nthetyp,nthetyp
259 do k=-nthetyp,nthetyp
260 read (ithep,'(6a)',end=111,err=111) res1
261 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
262 c VAR: aa0thet is variable describing the average value of Foureir
263 c VAR: expansion series
264 c VAR: aathet is foureir expansion in theta/2 angle for full formula
265 c VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
266 Condens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
267 read (ithep,*,end=111,err=111)
268 &(aathet(l,i,j,k,iblock),l=1,ntheterm)
269 read (ithep,*,end=111,err=111)
270 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
271 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
272 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
273 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),
275 read (ithep,*,end=111,err=111)
276 & (((ffthet(llll,lll,ll,i,j,k,iblock),
277 & ffthet(lll,llll,ll,i,j,k,iblock),
278 & ggthet(llll,lll,ll,i,j,k,iblock),
279 & ggthet(lll,llll,ll,i,j,k,iblock),
280 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
285 C For dummy ends assign glycine-type coefficients of theta-only terms; the
286 C coefficients of theta-and-gamma-dependent terms are zero.
287 C IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
288 C RECOMENTDED AFTER VERSION 3.3)
292 c aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
293 c aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
295 c aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
296 c aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
299 c aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
301 c aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
304 C AND COMMENT THE LOOPS BELOW
308 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
309 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
311 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
312 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
315 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
317 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
321 C Substitution for D aminoacids from symmetry.
324 do j=-nthetyp,nthetyp
325 do k=-nthetyp,nthetyp
326 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
328 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
332 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
333 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
334 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
335 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
341 ffthet(llll,lll,ll,i,j,k,iblock)=
342 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
343 ffthet(lll,llll,ll,i,j,k,iblock)=
344 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
345 ggthet(llll,lll,ll,i,j,k,iblock)=
346 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
347 ggthet(lll,llll,ll,i,j,k,iblock)=
348 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
357 C Control printout of the coefficients of virtual-bond-angle potentials
360 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
364 write (iout,'(//4a)')
365 & 'Type ',onelett(i),onelett(j),onelett(k)
366 write (iout,'(//a,10x,a)') " l","a[l]"
367 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
368 write (iout,'(i2,1pe15.5)')
369 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
371 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))')
372 & "b",l,"c",l,"d",l,"e",l
374 write (iout,'(i2,4(1pe15.5))') m,
375 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
376 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
380 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
381 & "f+",l,"f-",l,"g+",l,"g-",l
384 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
385 & ffthet(n,m,l,i,j,k,iblock),
386 & ffthet(m,n,l,i,j,k,iblock),
387 & ggthet(n,m,l,i,j,k,iblock),
388 & ggthet(m,n,l,i,j,k,iblock)
397 write (2,*) "Start reading THETA_PDB",ithep_pdb
400 read (ithep_pdb,*,err=111,end=111)
401 & a0thet(i),(athet(j,i,1,1),j=1,2),
402 & (bthet(j,i,1,1),j=1,2)
403 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
404 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
405 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
409 athet(1,i,1,-1)=athet(1,i,1,1)
410 athet(2,i,1,-1)=athet(2,i,1,1)
411 bthet(1,i,1,-1)=-bthet(1,i,1,1)
412 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)
420 athet(1,i,-1,-1)=athet(1,-i,1,1)
421 athet(2,i,-1,-1)=-athet(2,-i,1,1)
422 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
423 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
424 athet(1,i,-1,1)=athet(1,-i,1,1)
425 athet(2,i,-1,1)=-athet(2,-i,1,1)
426 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
427 bthet(2,i,-1,1)=bthet(2,-i,1,1)
428 athet(1,i,1,-1)=-athet(1,-i,1,1)
429 athet(2,i,1,-1)=athet(2,-i,1,1)
430 bthet(1,i,1,-1)=bthet(1,-i,1,1)
431 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
436 polthet(j,i)=polthet(j,-i)
439 gthet(j,i)=gthet(j,-i)
442 write (2,*) "End reading THETA_PDB"
448 C Read the parameters of the probability distribution/energy expression
449 C of the side chains.
452 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
456 dsc_inv(i)=1.0D0/dsc(i)
467 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
468 & ((blower(k,l,1),l=1,k),k=1,3)
469 censc(1,1,-i)=censc(1,1,i)
470 censc(2,1,-i)=censc(2,1,i)
471 censc(3,1,-i)=-censc(3,1,i)
473 read (irotam,*,end=112,err=112) bsc(j,i)
474 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
475 & ((blower(k,l,j),l=1,k),k=1,3)
476 censc(1,j,-i)=censc(1,j,i)
477 censc(2,j,-i)=censc(2,j,i)
478 censc(3,j,-i)=-censc(3,j,i)
479 C BSC is amplitude of Gaussian
486 akl=akl+blower(k,m,j)*blower(l,m,j)
490 if (((k.eq.3).and.(l.ne.3))
491 & .or.((l.eq.3).and.(k.ne.3))) then
492 gaussc(k,l,j,-i)=-akl
493 gaussc(l,k,j,-i)=-akl
505 write (iout,'(/a)') 'Parameters of side-chain local geometry'
510 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
511 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
512 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
513 & 'log h',(bsc(j,i),j=1,nlobi)
514 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
515 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
517 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
518 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
521 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
522 write (iout,'(a,f10.4,4(16x,f10.4))')
523 & 'Center ',(bsc(j,i),j=1,nlobi)
524 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
533 C Read scrot parameters for potentials determined from all-atom AM1 calculations
534 C added by Urszula Kozlowska 07/11/2007
537 read (irotam,*,end=112,err=112)
539 read (irotam,*,end=112,err=112)
542 read(irotam,*,end=112,err=112) sc_parmin(j,i)
547 C Read the parameters of the probability distribution/energy expression
548 C of the side chains.
550 write (2,*) "Start reading ROTAM_PDB"
552 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
556 dsc_inv(i)=1.0D0/dsc(i)
567 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
568 & ((blower(k,l,1),l=1,k),k=1,3)
570 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
571 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
572 & ((blower(k,l,j),l=1,k),k=1,3)
579 akl=akl+blower(k,m,j)*blower(l,m,j)
589 write (2,*) "End reading ROTAM_PDB"
595 C Read torsional parameters in old format
597 read (itorp,*,end=113,err=113) ntortyp,nterm_old
598 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
599 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
604 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
610 write (iout,'(/a/)') 'Torsional constants:'
613 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
614 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
620 C Read torsional parameters
622 read (itorp,*,end=113,err=113) ntortyp
623 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
626 itortyp(i)=-itortyp(-i)
628 write (iout,*) 'ntortyp',ntortyp
630 do j=-ntortyp+1,ntortyp-1
631 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
633 nterm(-i,-j,iblock)=nterm(i,j,iblock)
634 nlor(-i,-j,iblock)=nlor(i,j,iblock)
637 do k=1,nterm(i,j,iblock)
638 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
640 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
641 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
642 v0ij=v0ij+si*v1(k,i,j,iblock)
644 c write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
645 c write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
646 c &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
648 do k=1,nlor(i,j,iblock)
649 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
650 & vlor2(k,i,j),vlor3(k,i,j)
651 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
654 v0(-i,-j,iblock)=v0ij
660 write (iout,'(/a/)') 'Torsional constants:'
663 write (iout,*) 'ityp',i,' jtyp',j
664 write (iout,*) 'Fourier constants'
665 do k=1,nterm(i,j,iblock)
666 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
669 write (iout,*) 'Lorenz constants'
670 do k=1,nlor(i,j,iblock)
671 write (iout,'(3(1pe15.5))')
672 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
679 C 6/23/01 Read parameters for double torsionals
683 do j=-ntortyp+1,ntortyp-1
684 do k=-ntortyp+1,ntortyp-1
685 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
686 c write (iout,*) "OK onelett",
689 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
690 & .or. t3.ne.toronelet(k)) then
691 write (iout,*) "Error in double torsional parameter file",
694 call MPI_Finalize(Ierror)
696 stop "Error in double torsional parameter file"
698 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
699 & ntermd_2(i,j,k,iblock)
700 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
701 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
702 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
703 & ntermd_1(i,j,k,iblock))
704 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
705 & ntermd_1(i,j,k,iblock))
706 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
707 & ntermd_1(i,j,k,iblock))
708 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
709 & ntermd_1(i,j,k,iblock))
710 C Martix of D parameters for one dimesional foureir series
711 do l=1,ntermd_1(i,j,k,iblock)
712 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
713 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
714 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
715 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
716 c write(iout,*) "whcodze" ,
717 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
719 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
720 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
721 & v2s(m,l,i,j,k,iblock),
722 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
723 C Martix of D parameters for two dimesional fourier series
724 do l=1,ntermd_2(i,j,k,iblock)
726 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
727 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
728 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
729 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
738 write (iout,*) 'Constants for double torsionals'
741 do j=-ntortyp+1,ntortyp-1
742 do k=-ntortyp+1,ntortyp-1
743 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
744 & ' nsingle',ntermd_1(i,j,k,iblock),
745 & ' ndouble',ntermd_2(i,j,k,iblock)
747 write (iout,*) 'Single angles:'
748 do l=1,ntermd_1(i,j,k,iblock)
749 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
750 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
751 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
752 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
755 write (iout,*) 'Pairs of angles:'
756 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
757 do l=1,ntermd_2(i,j,k,iblock)
758 write (iout,'(i5,20f10.5)')
759 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
762 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
763 do l=1,ntermd_2(i,j,k,iblock)
764 write (iout,'(i5,20f10.5)')
765 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
766 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
775 C Read of Side-chain backbone correlation parameters
776 C Modified 11 May 2012 by Adasko
779 read (isccor,*,end=119,err=119) nsccortyp
781 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
783 isccortyp(i)=-isccortyp(-i)
785 iscprol=isccortyp(20)
786 c write (iout,*) 'ntortyp',ntortyp
788 cc maxinter is maximum interaction sites
792 read (isccor,*,end=119,err=119)
793 &nterm_sccor(i,j),nlor_sccor(i,j)
799 nterm_sccor(-i,j)=nterm_sccor(i,j)
800 nterm_sccor(-i,-j)=nterm_sccor(i,j)
801 nterm_sccor(i,-j)=nterm_sccor(i,j)
802 do k=1,nterm_sccor(i,j)
803 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
805 if (j.eq.iscprol) then
806 if (i.eq.isccortyp(10)) then
807 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
808 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
810 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
811 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
812 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
813 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
814 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
815 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
816 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
817 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
820 if (i.eq.isccortyp(10)) then
821 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
822 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
824 if (j.eq.isccortyp(10)) then
825 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
826 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
828 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
829 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
830 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
831 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
832 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
833 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
837 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
838 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
839 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
840 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
843 do k=1,nlor_sccor(i,j)
844 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
845 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
846 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
847 &(1+vlor3sccor(k,i,j)**2)
849 v0sccor(l,i,j)=v0ijsccor
850 v0sccor(l,-i,j)=v0ijsccor1
851 v0sccor(l,i,-j)=v0ijsccor2
852 v0sccor(l,-i,-j)=v0ijsccor3
858 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
859 c write (iout,*) 'ntortyp',ntortyp
861 cc maxinter is maximum interaction sites
865 read (isccor,*,end=119,err=119)
866 & nterm_sccor(i,j),nlor_sccor(i,j)
870 do k=1,nterm_sccor(i,j)
871 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
873 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
876 do k=1,nlor_sccor(i,j)
877 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
878 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
879 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
880 &(1+vlor3sccor(k,i,j)**2)
882 v0sccor(i,j,iblock)=v0ijsccor
890 write (iout,'(/a/)') 'Torsional constants:'
893 write (iout,*) 'ityp',i,' jtyp',j
894 write (iout,*) 'Fourier constants'
895 do k=1,nterm_sccor(i,j)
896 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
898 write (iout,*) 'Lorenz constants'
899 do k=1,nlor_sccor(i,j)
900 write (iout,'(3(1pe15.5))')
901 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
908 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
909 C interaction energy of the Gly, Ala, and Pro prototypes.
913 write (iout,*) "Coefficients of the cumulants"
915 read (ifourier,*) nloctyp
917 read (ifourier,*,end=115,err=115)
918 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
920 write (iout,*) 'Type',i
921 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
963 c Ctilde(1,1,i)=0.0d0
964 c Ctilde(1,2,i)=0.0d0
965 c Ctilde(2,1,i)=0.0d0
966 c Ctilde(2,2,i)=0.0d0
988 c Dtilde(1,1,i)=0.0d0
989 c Dtilde(1,2,i)=0.0d0
990 c Dtilde(2,1,i)=0.0d0
991 c Dtilde(2,2,i)=0.0d0
992 EE(1,1,i)= b(10)+b(11)
993 EE(2,2,i)=-b(10)+b(11)
994 EE(2,1,i)= b(12)-b(13)
995 EE(1,2,i)= b(12)+b(13)
996 EE(1,1,-i)= b(10)+b(11)
997 EE(2,2,-i)=-b(10)+b(11)
998 EE(2,1,-i)=-b(12)+b(13)
999 EE(1,2,-i)=-b(12)-b(13)
1005 c ee(2,1,i)=ee(1,2,i)
1009 write (iout,*) 'Type',i
1011 write(iout,*) B1(1,i),B1(2,i)
1013 write(iout,*) B2(1,i),B2(2,i)
1016 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1020 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1024 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1030 C Read electrostatic-interaction parameters
1034 write (iout,'(/a)') 'Electrostatic interaction constants:'
1035 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1036 & 'IT','JT','APP','BPP','AEL6','AEL3'
1038 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1039 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1040 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1041 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1046 app (i,j)=epp(i,j)*rri*rri
1047 bpp (i,j)=-2.0D0*epp(i,j)*rri
1048 ael6(i,j)=elpp6(i,j)*4.2D0**6
1049 ael3(i,j)=elpp3(i,j)*4.2D0**3
1051 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1052 & ael6(i,j),ael3(i,j)
1057 C Read side-chain interaction parameters.
1059 read (isidep,*,end=117,err=117) ipot,expon
1060 if (ipot.lt.1 .or. ipot.gt.5) then
1061 write (iout,'(2a)') 'Error while reading SC interaction',
1062 & 'potential file - unknown potential type.'
1064 call MPI_Finalize(Ierror)
1070 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1071 & ', exponents are ',expon,2*expon
1072 goto (10,20,30,30,40) ipot
1073 C----------------------- LJ potential ---------------------------------
1074 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1075 & (sigma0(i),i=1,ntyp)
1077 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1078 write (iout,'(a/)') 'The epsilon array:'
1079 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1080 write (iout,'(/a)') 'One-body parameters:'
1081 write (iout,'(a,4x,a)') 'residue','sigma'
1082 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1085 C----------------------- LJK potential --------------------------------
1086 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1087 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1089 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1090 write (iout,'(a/)') 'The epsilon array:'
1091 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1092 write (iout,'(/a)') 'One-body parameters:'
1093 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1094 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1098 C---------------------- GB or BP potential -----------------------------
1100 read (isidep,*,end=116,err=116)(eps(i,j),j=i,ntyp)
1102 read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
1103 read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
1104 read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
1105 read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
1106 C For the GB potential convert sigma'**2 into chi'
1109 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1113 write (iout,'(/a/)') 'Parameters of the BP potential:'
1114 write (iout,'(a/)') 'The epsilon array:'
1115 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1116 write (iout,'(/a)') 'One-body parameters:'
1117 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1119 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1120 & chip(i),alp(i),i=1,ntyp)
1123 C--------------------- GBV potential -----------------------------------
1124 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1125 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1126 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1128 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1129 write (iout,'(a/)') 'The epsilon array:'
1130 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1131 write (iout,'(/a)') 'One-body parameters:'
1132 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1133 & 's||/s_|_^2',' chip ',' alph '
1134 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1135 & sigii(i),chip(i),alp(i),i=1,ntyp)
1139 C-----------------------------------------------------------------------
1140 C Calculate the "working" parameters of SC interactions.
1148 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1149 sigma(j,i)=sigma(i,j)
1150 rs0(i,j)=dwa16*sigma(i,j)
1154 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1155 & 'Working parameters of the SC interactions:',
1156 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1161 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1170 sigeps=dsign(1.0D0,epsij)
1172 aa(i,j)=epsij*rrij*rrij
1173 bb(i,j)=-sigeps*epsij*rrij
1177 sigt1sq=sigma0(i)**2
1178 sigt2sq=sigma0(j)**2
1181 ratsig1=sigt2sq/sigt1sq
1182 ratsig2=1.0D0/ratsig1
1183 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1184 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1185 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1189 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1190 sigmaii(i,j)=rsum_max
1191 sigmaii(j,i)=rsum_max
1193 c sigmaii(i,j)=r0(i,j)
1194 c sigmaii(j,i)=r0(i,j)
1196 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1197 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1198 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1199 augm(i,j)=epsij*r_augm**(2*expon)
1200 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1207 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1208 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1209 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1215 C Define the SC-p interaction constants (hard-coded; old style)
1218 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1220 c aad(i,1)=0.3D0*4.0D0**12
1221 C Following line for constants currently implemented
1222 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1223 aad(i,1)=1.5D0*4.0D0**12
1224 c aad(i,1)=0.17D0*5.6D0**12
1226 C "Soft" SC-p repulsion
1228 C Following line for constants currently implemented
1229 c aad(i,1)=0.3D0*4.0D0**6
1230 C "Hard" SC-p repulsion
1231 bad(i,1)=3.0D0*4.0D0**6
1232 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1241 C 8/9/01 Read the SC-p interaction constants from file
1244 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1247 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1248 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1249 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1250 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1254 write (iout,*) "Parameters of SC-p interactions:"
1256 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1257 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1263 C Define the constants of the disulfide bridge
1267 c Old arbitrary potential - commented out.
1272 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1273 c energy surface of diethyl disulfide.
1274 c A. Liwo and U. Kozlowska, 11/24/03
1291 write (iout,'(/a)') "Disulfide bridge parameters:"
1292 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1293 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1294 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1295 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1299 111 write (iout,*) "Error reading bending energy parameters."
1301 112 write (iout,*) "Error reading rotamer energy parameters."
1303 113 write (iout,*) "Error reading torsional energy parameters."
1305 114 write (iout,*) "Error reading double torsional energy parameters."
1308 & "Error reading cumulant (multibody energy) parameters."
1310 116 write (iout,*) "Error reading electrostatic energy parameters."
1312 117 write (iout,*) "Error reading side chain interaction parameters."
1314 118 write (iout,*) "Error reading SCp interaction parameters."
1316 119 write (iout,*) "Error reading SCCOR parameters"
1319 call MPI_Finalize(Ierror)
1326 subroutine getenv_loc(var, val)
1327 character(*) var, val
1330 character(2000) line
1333 open (196,file='env',status='old',readonly,shared)
1335 c write(*,*)'looking for ',var
1336 10 read(196,*,err=11,end=11)line
1337 iread=index(line,var)
1338 c write(*,*)iread,' ',var,' ',line
1339 if (iread.eq.0) go to 10
1340 c write(*,*)'---> ',line
1346 iread=iread+ilen(var)+1
1347 read (line(iread:),*,err=12,end=12) val
1348 c write(*,*)'OK: ',var,' = ',val
1354 #elif (defined CRAY)
1355 integer lennam,lenval,ierror
1357 c getenv using a POSIX call, useful on the T3D
1358 c Sept 1996, comment out error check on advice of H. Pritchard
1361 if(lennam.le.0) stop '--error calling getenv--'
1362 call pxfgetenv(var,lennam,val,lenval,ierror)
1363 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1365 call getenv(var,val)