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 c 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 c 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 c 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 c 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 c 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
918 read (ifourier,*,end=115,err=115)
919 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
921 read (ifourier,*,end=115,err=115) (bnew1(ii,1,i),ii=1,3)
922 read (ifourier,*,end=115,err=115) (bnew2(ii,1,i),ii=1,3)
923 read (ifourier,*,end=115,err=115) (bnew1(ii,2,i),ii=1,1)
924 read (ifourier,*,end=115,err=115) (bnew2(ii,2,i),ii=1,1)
925 read (ifourier,*,end=115,err=115) (eenew(ii,i),ii=1,1)
928 write (iout,*) 'Type',i
929 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
937 c B1tilde(1,i) = b(3)
938 c B1tilde(2,i) =-b(5)
939 c B1tilde(1,-i) =-b(3)
940 c B1tilde(2,-i) =b(5)
964 Ctilde(2,1,i)=-b(9,i)
966 Ctilde(1,1,-i)=b(7,i)
967 Ctilde(1,2,-i)=-b(9,i)
968 Ctilde(2,1,-i)=b(9,i)
969 Ctilde(2,2,-i)=b(7,i)
971 c Ctilde(1,1,i)=0.0d0
972 c Ctilde(1,2,i)=0.0d0
973 c Ctilde(2,1,i)=0.0d0
974 c Ctilde(2,2,i)=0.0d0
989 Dtilde(2,1,i)=-b(8,i)
991 Dtilde(1,1,-i)=b(6,i)
992 Dtilde(1,2,-i)=-b(8,i)
993 Dtilde(2,1,-i)=b(8,i)
994 Dtilde(2,2,-i)=b(6,i)
996 c Dtilde(1,1,i)=0.0d0
997 c Dtilde(1,2,i)=0.0d0
998 c Dtilde(2,1,i)=0.0d0
999 c Dtilde(2,2,i)=0.0d0
1000 EEold(1,1,i)= b(10,i)+b(11,i)
1001 EEold(2,2,i)=-b(10,i)+b(11,i)
1002 EEold(2,1,i)= b(12,i)-b(13,i)
1003 EEold(1,2,i)= b(12,i)+b(13,i)
1004 EEold(1,1,-i)= b(10,i)+b(11,i)
1005 EEold(2,2,-i)=-b(10,i)+b(11,i)
1006 EEold(2,1,-i)=-b(12,i)+b(13,i)
1007 EEold(1,2,-i)=-b(12,i)-b(13,i)
1012 c ee(2,1,i)=ee(1,2,i)
1017 write (iout,*) 'Type',i
1019 write(iout,*) B1(1,i),B1(2,i)
1021 write(iout,*) B2(1,i),B2(2,i)
1024 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1028 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1032 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1039 C Read electrostatic-interaction parameters
1043 write (iout,'(/a)') 'Electrostatic interaction constants:'
1044 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1045 & 'IT','JT','APP','BPP','AEL6','AEL3'
1047 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1048 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1049 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1050 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1055 app (i,j)=epp(i,j)*rri*rri
1056 bpp (i,j)=-2.0D0*epp(i,j)*rri
1057 ael6(i,j)=elpp6(i,j)*4.2D0**6
1058 ael3(i,j)=elpp3(i,j)*4.2D0**3
1060 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1061 & ael6(i,j),ael3(i,j)
1066 C Read side-chain interaction parameters.
1068 read (isidep,*,end=117,err=117) ipot,expon
1069 if (ipot.lt.1 .or. ipot.gt.5) then
1070 write (iout,'(2a)') 'Error while reading SC interaction',
1071 & 'potential file - unknown potential type.'
1073 call MPI_Finalize(Ierror)
1079 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1080 & ', exponents are ',expon,2*expon
1081 goto (10,20,30,30,40) ipot
1082 C----------------------- LJ potential ---------------------------------
1083 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1084 & (sigma0(i),i=1,ntyp)
1086 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1087 write (iout,'(a/)') 'The epsilon array:'
1088 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1089 write (iout,'(/a)') 'One-body parameters:'
1090 write (iout,'(a,4x,a)') 'residue','sigma'
1091 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1094 C----------------------- LJK potential --------------------------------
1095 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1096 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1098 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1099 write (iout,'(a/)') 'The epsilon array:'
1100 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1101 write (iout,'(/a)') 'One-body parameters:'
1102 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1103 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1107 C---------------------- GB or BP potential -----------------------------
1109 read (isidep,*,end=116,err=116)(eps(i,j),j=i,ntyp)
1111 read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
1112 read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
1113 read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
1114 read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
1115 C now we start reading lipid
1117 read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
1118 C print *,"WARNING!!"
1120 C epslip(i,j)=epslip(i,j)+0.05d0
1123 c write(iout,*) epslip(1,1),"OK?"
1124 C For the GB potential convert sigma'**2 into chi'
1127 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1131 write (iout,'(/a/)') 'Parameters of the BP potential:'
1132 write (iout,'(a/)') 'The epsilon array:'
1133 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1134 write (iout,'(/a)') 'One-body parameters:'
1135 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1137 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1138 & chip(i),alp(i),i=1,ntyp)
1141 C--------------------- GBV potential -----------------------------------
1142 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1143 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1144 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1146 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1147 write (iout,'(a/)') 'The epsilon array:'
1148 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1149 write (iout,'(/a)') 'One-body parameters:'
1150 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1151 & 's||/s_|_^2',' chip ',' alph '
1152 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1153 & sigii(i),chip(i),alp(i),i=1,ntyp)
1157 C-----------------------------------------------------------------------
1158 C Calculate the "working" parameters of SC interactions.
1162 epslip(i,j)=epslip(j,i)
1167 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1168 sigma(j,i)=sigma(i,j)
1169 rs0(i,j)=dwa16*sigma(i,j)
1173 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1174 & 'Working parameters of the SC interactions:',
1175 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1180 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1189 sigeps=dsign(1.0D0,epsij)
1191 aa_aq(i,j)=epsij*rrij*rrij
1192 bb_aq(i,j)=-sigeps*epsij*rrij
1193 aa_aq(j,i)=aa_aq(i,j)
1194 bb_aq(j,i)=bb_aq(i,j)
1195 epsijlip=epslip(i,j)
1196 sigeps=dsign(1.0D0,epsijlip)
1197 epsijlip=dabs(epsijlip)
1198 aa_lip(i,j)=epsijlip*rrij*rrij
1199 bb_lip(i,j)=-sigeps*epsijlip*rrij
1200 aa_lip(j,i)=aa_lip(i,j)
1201 bb_lip(j,i)=bb_lip(i,j)
1203 sigt1sq=sigma0(i)**2
1204 sigt2sq=sigma0(j)**2
1207 ratsig1=sigt2sq/sigt1sq
1208 ratsig2=1.0D0/ratsig1
1209 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1210 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1211 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1215 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1216 sigmaii(i,j)=rsum_max
1217 sigmaii(j,i)=rsum_max
1219 c sigmaii(i,j)=r0(i,j)
1220 c sigmaii(j,i)=r0(i,j)
1222 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1223 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1224 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1225 augm(i,j)=epsij*r_augm**(2*expon)
1226 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1233 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1234 & restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1235 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1241 C Define the SC-p interaction constants (hard-coded; old style)
1244 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1246 c aad(i,1)=0.3D0*4.0D0**12
1247 C Following line for constants currently implemented
1248 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1249 aad(i,1)=1.5D0*4.0D0**12
1250 c aad(i,1)=0.17D0*5.6D0**12
1252 C "Soft" SC-p repulsion
1254 C Following line for constants currently implemented
1255 c aad(i,1)=0.3D0*4.0D0**6
1256 C "Hard" SC-p repulsion
1257 bad(i,1)=3.0D0*4.0D0**6
1258 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1267 C 8/9/01 Read the SC-p interaction constants from file
1270 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1273 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1274 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1275 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1276 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1280 write (iout,*) "Parameters of SC-p interactions:"
1282 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1283 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1289 C Define the constants of the disulfide bridge
1293 c Old arbitrary potential - commented out.
1298 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1299 c energy surface of diethyl disulfide.
1300 c A. Liwo and U. Kozlowska, 11/24/03
1317 write (iout,'(/a)') "Disulfide bridge parameters:"
1318 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1319 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1320 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1321 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1325 111 write (iout,*) "Error reading bending energy parameters."
1327 112 write (iout,*) "Error reading rotamer energy parameters."
1329 113 write (iout,*) "Error reading torsional energy parameters."
1331 114 write (iout,*) "Error reading double torsional energy parameters."
1334 & "Error reading cumulant (multibody energy) parameters."
1336 116 write (iout,*) "Error reading electrostatic energy parameters."
1338 1161 write (iout,*) "Error reading electrostatic energy parameters.Lip"
1340 117 write (iout,*) "Error reading side chain interaction parameters."
1342 118 write (iout,*) "Error reading SCp interaction parameters."
1344 119 write (iout,*) "Error reading SCCOR parameters"
1347 call MPI_Finalize(Ierror)
1354 subroutine getenv_loc(var, val)
1355 character(*) var, val
1358 character(2000) line
1361 open (196,file='env',status='old',readonly,shared)
1363 c write(*,*)'looking for ',var
1364 10 read(196,*,err=11,end=11)line
1365 iread=index(line,var)
1366 c write(*,*)iread,' ',var,' ',line
1367 if (iread.eq.0) go to 10
1368 c write(*,*)'---> ',line
1374 iread=iread+ilen(var)+1
1375 read (line(iread:),*,err=12,end=12) val
1376 c write(*,*)'OK: ',var,' = ',val
1382 #elif (defined CRAY)
1383 integer lennam,lenval,ierror
1385 c getenv using a POSIX call, useful on the T3D
1386 c Sept 1996, comment out error check on advice of H. Pritchard
1389 if(lennam.le.0) stop '--error calling getenv--'
1390 call pxfgetenv(var,lennam,val,lenval,ierror)
1391 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1393 call getenv(var,val)