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 now we start reading lipid
1108 read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
1111 epslip(i,j)=epslip(i,j)+0.5d0
1114 C For the GB potential convert sigma'**2 into chi'
1117 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1121 write (iout,'(/a/)') 'Parameters of the BP potential:'
1122 write (iout,'(a/)') 'The epsilon array:'
1123 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1124 write (iout,'(/a)') 'One-body parameters:'
1125 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1127 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1128 & chip(i),alp(i),i=1,ntyp)
1131 C--------------------- GBV potential -----------------------------------
1132 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1133 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1134 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1136 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1137 write (iout,'(a/)') 'The epsilon array:'
1138 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1139 write (iout,'(/a)') 'One-body parameters:'
1140 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1141 & 's||/s_|_^2',' chip ',' alph '
1142 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1143 & sigii(i),chip(i),alp(i),i=1,ntyp)
1147 C-----------------------------------------------------------------------
1148 C Calculate the "working" parameters of SC interactions.
1156 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1157 sigma(j,i)=sigma(i,j)
1158 rs0(i,j)=dwa16*sigma(i,j)
1162 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1163 & 'Working parameters of the SC interactions:',
1164 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1169 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1178 sigeps=dsign(1.0D0,epsij)
1180 aa_aq(i,j)=epsij*rrij*rrij
1181 bb_aq(i,j)=-sigeps*epsij*rrij
1182 aa_aq(j,i)=aa_aq(i,j)
1183 bb_aq(j,i)=bb_aq(i,j)
1184 epsijlip=epslip(i,j)
1185 sigeps=dsign(1.0D0,epsijlip)
1186 epsijlip=dabs(epsijlip)
1187 aa_lip(i,j)=epsijlip*rrij*rrij
1188 bb_lip(i,j)=-sigeps*epsijlip*rrij
1189 aa_lip(j,i)=aa_lip(i,j)
1190 bb_lip(j,i)=bb_lip(i,j)
1192 sigt1sq=sigma0(i)**2
1193 sigt2sq=sigma0(j)**2
1196 ratsig1=sigt2sq/sigt1sq
1197 ratsig2=1.0D0/ratsig1
1198 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1199 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1200 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1204 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1205 sigmaii(i,j)=rsum_max
1206 sigmaii(j,i)=rsum_max
1208 c sigmaii(i,j)=r0(i,j)
1209 c sigmaii(j,i)=r0(i,j)
1211 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1212 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1213 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1214 augm(i,j)=epsij*r_augm**(2*expon)
1215 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1222 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1223 & restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1224 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1230 C Define the SC-p interaction constants (hard-coded; old style)
1233 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1235 c aad(i,1)=0.3D0*4.0D0**12
1236 C Following line for constants currently implemented
1237 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1238 aad(i,1)=1.5D0*4.0D0**12
1239 c aad(i,1)=0.17D0*5.6D0**12
1241 C "Soft" SC-p repulsion
1243 C Following line for constants currently implemented
1244 c aad(i,1)=0.3D0*4.0D0**6
1245 C "Hard" SC-p repulsion
1246 bad(i,1)=3.0D0*4.0D0**6
1247 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1256 C 8/9/01 Read the SC-p interaction constants from file
1259 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1262 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1263 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1264 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1265 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1269 write (iout,*) "Parameters of SC-p interactions:"
1271 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1272 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1278 C Define the constants of the disulfide bridge
1282 c Old arbitrary potential - commented out.
1287 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1288 c energy surface of diethyl disulfide.
1289 c A. Liwo and U. Kozlowska, 11/24/03
1306 write (iout,'(/a)') "Disulfide bridge parameters:"
1307 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1308 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1309 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1310 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1314 111 write (iout,*) "Error reading bending energy parameters."
1316 112 write (iout,*) "Error reading rotamer energy parameters."
1318 113 write (iout,*) "Error reading torsional energy parameters."
1320 114 write (iout,*) "Error reading double torsional energy parameters."
1323 & "Error reading cumulant (multibody energy) parameters."
1325 116 write (iout,*) "Error reading electrostatic energy parameters."
1327 1161 write (iout,*) "Error reading electrostatic energy parameters.Lip"
1329 117 write (iout,*) "Error reading side chain interaction parameters."
1331 118 write (iout,*) "Error reading SCp interaction parameters."
1333 119 write (iout,*) "Error reading SCCOR parameters"
1336 call MPI_Finalize(Ierror)
1343 subroutine getenv_loc(var, val)
1344 character(*) var, val
1347 character(2000) line
1350 open (196,file='env',status='old',readonly,shared)
1352 c write(*,*)'looking for ',var
1353 10 read(196,*,err=11,end=11)line
1354 iread=index(line,var)
1355 c write(*,*)iread,' ',var,' ',line
1356 if (iread.eq.0) go to 10
1357 c write(*,*)'---> ',line
1363 iread=iread+ilen(var)+1
1364 read (line(iread:),*,err=12,end=12) val
1365 c write(*,*)'OK: ',var,' = ',val
1371 #elif (defined CRAY)
1372 integer lennam,lenval,ierror
1374 c getenv using a POSIX call, useful on the T3D
1375 c Sept 1996, comment out error check on advice of H. Pritchard
1378 if(lennam.le.0) stop '--error calling getenv--'
1379 call pxfgetenv(var,lennam,val,lenval,ierror)
1380 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1382 call getenv(var,val)