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)
941 B1tilde(1,i) = b(3,i)
942 B1tilde(2,i) =-b(5,i)
943 B1tilde(1,-i) =-b(3,i)
944 B1tilde(2,-i) =b(5,i)
967 Ctilde(2,1,i)=-b(9,i)
969 Ctilde(1,1,-i)=b(7,i)
970 Ctilde(1,2,-i)=-b(9,i)
971 Ctilde(2,1,-i)=b(9,i)
972 Ctilde(2,2,-i)=b(7,i)
974 c Ctilde(1,1,i)=0.0d0
975 c Ctilde(1,2,i)=0.0d0
976 c Ctilde(2,1,i)=0.0d0
977 c Ctilde(2,2,i)=0.0d0
992 Dtilde(2,1,i)=-b(8,i)
994 Dtilde(1,1,-i)=b(6,i)
995 Dtilde(1,2,-i)=-b(8,i)
996 Dtilde(2,1,-i)=b(8,i)
997 Dtilde(2,2,-i)=b(6,i)
999 c Dtilde(1,1,i)=0.0d0
1000 c Dtilde(1,2,i)=0.0d0
1001 c Dtilde(2,1,i)=0.0d0
1002 c Dtilde(2,2,i)=0.0d0
1003 EEold(1,1,i)= b(10,i)+b(11,i)
1004 EEold(2,2,i)=-b(10,i)+b(11,i)
1005 EEold(2,1,i)= b(12,i)-b(13,i)
1006 EEold(1,2,i)= b(12,i)+b(13,i)
1007 EEold(1,1,-i)= b(10,i)+b(11,i)
1008 EEold(2,2,-i)=-b(10,i)+b(11,i)
1009 EEold(2,1,-i)=-b(12,i)+b(13,i)
1010 EEold(1,2,-i)=-b(12,i)-b(13,i)
1015 c ee(2,1,i)=ee(1,2,i)
1020 write (iout,*) 'Type',i
1022 write(iout,*) B1(1,i),B1(2,i)
1024 write(iout,*) B2(1,i),B2(2,i)
1027 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1031 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1035 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1042 C Read electrostatic-interaction parameters
1046 write (iout,'(/a)') 'Electrostatic interaction constants:'
1047 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1048 & 'IT','JT','APP','BPP','AEL6','AEL3'
1050 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1051 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1052 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1053 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1058 app (i,j)=epp(i,j)*rri*rri
1059 bpp (i,j)=-2.0D0*epp(i,j)*rri
1060 ael6(i,j)=elpp6(i,j)*4.2D0**6
1061 ael3(i,j)=elpp3(i,j)*4.2D0**3
1063 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1064 & ael6(i,j),ael3(i,j)
1069 C Read side-chain interaction parameters.
1071 read (isidep,*,end=117,err=117) ipot,expon
1072 if (ipot.lt.1 .or. ipot.gt.5) then
1073 write (iout,'(2a)') 'Error while reading SC interaction',
1074 & 'potential file - unknown potential type.'
1076 call MPI_Finalize(Ierror)
1082 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1083 & ', exponents are ',expon,2*expon
1084 goto (10,20,30,30,40) ipot
1085 C----------------------- LJ potential ---------------------------------
1086 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1087 & (sigma0(i),i=1,ntyp)
1089 write (iout,'(/a/)') 'Parameters of the LJ 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,a)') 'residue','sigma'
1094 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1097 C----------------------- LJK potential --------------------------------
1098 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1099 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1101 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1102 write (iout,'(a/)') 'The epsilon array:'
1103 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1104 write (iout,'(/a)') 'One-body parameters:'
1105 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1106 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1110 C---------------------- GB or BP potential -----------------------------
1112 read (isidep,*,end=116,err=116)(eps(i,j),j=i,ntyp)
1114 read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
1115 read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
1116 read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
1117 read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
1118 C now we start reading lipid
1120 read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
1121 C print *,"WARNING!!"
1123 C epslip(i,j)=epslip(i,j)+0.05d0
1126 c write(iout,*) epslip(1,1),"OK?"
1127 C For the GB potential convert sigma'**2 into chi'
1130 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1134 write (iout,'(/a/)') 'Parameters of the BP potential:'
1135 write (iout,'(a/)') 'The epsilon array:'
1136 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1137 write (iout,'(/a)') 'One-body parameters:'
1138 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1140 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1141 & chip(i),alp(i),i=1,ntyp)
1144 C--------------------- GBV potential -----------------------------------
1145 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1146 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1147 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1149 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1150 write (iout,'(a/)') 'The epsilon array:'
1151 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1152 write (iout,'(/a)') 'One-body parameters:'
1153 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1154 & 's||/s_|_^2',' chip ',' alph '
1155 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1156 & sigii(i),chip(i),alp(i),i=1,ntyp)
1160 C-----------------------------------------------------------------------
1161 C Calculate the "working" parameters of SC interactions.
1165 epslip(i,j)=epslip(j,i)
1170 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1171 sigma(j,i)=sigma(i,j)
1172 rs0(i,j)=dwa16*sigma(i,j)
1176 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1177 & 'Working parameters of the SC interactions:',
1178 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1183 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1192 sigeps=dsign(1.0D0,epsij)
1194 aa_aq(i,j)=epsij*rrij*rrij
1195 bb_aq(i,j)=-sigeps*epsij*rrij
1196 aa_aq(j,i)=aa_aq(i,j)
1197 bb_aq(j,i)=bb_aq(i,j)
1198 epsijlip=epslip(i,j)
1199 sigeps=dsign(1.0D0,epsijlip)
1200 epsijlip=dabs(epsijlip)
1201 aa_lip(i,j)=epsijlip*rrij*rrij
1202 bb_lip(i,j)=-sigeps*epsijlip*rrij
1203 aa_lip(j,i)=aa_lip(i,j)
1204 bb_lip(j,i)=bb_lip(i,j)
1206 sigt1sq=sigma0(i)**2
1207 sigt2sq=sigma0(j)**2
1210 ratsig1=sigt2sq/sigt1sq
1211 ratsig2=1.0D0/ratsig1
1212 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1213 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1214 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1218 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1219 sigmaii(i,j)=rsum_max
1220 sigmaii(j,i)=rsum_max
1222 c sigmaii(i,j)=r0(i,j)
1223 c sigmaii(j,i)=r0(i,j)
1225 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1226 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1227 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1228 augm(i,j)=epsij*r_augm**(2*expon)
1229 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1236 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1237 & restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1238 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1244 C Define the SC-p interaction constants (hard-coded; old style)
1247 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1249 c aad(i,1)=0.3D0*4.0D0**12
1250 C Following line for constants currently implemented
1251 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1252 aad(i,1)=1.5D0*4.0D0**12
1253 c aad(i,1)=0.17D0*5.6D0**12
1255 C "Soft" SC-p repulsion
1257 C Following line for constants currently implemented
1258 c aad(i,1)=0.3D0*4.0D0**6
1259 C "Hard" SC-p repulsion
1260 bad(i,1)=3.0D0*4.0D0**6
1261 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1270 C 8/9/01 Read the SC-p interaction constants from file
1273 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1276 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1277 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1278 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1279 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1283 write (iout,*) "Parameters of SC-p interactions:"
1285 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1286 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1292 C Define the constants of the disulfide bridge
1296 c Old arbitrary potential - commented out.
1301 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1302 c energy surface of diethyl disulfide.
1303 c A. Liwo and U. Kozlowska, 11/24/03
1320 write (iout,'(/a)') "Disulfide bridge parameters:"
1321 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1322 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1323 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1324 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1328 111 write (iout,*) "Error reading bending energy parameters."
1330 112 write (iout,*) "Error reading rotamer energy parameters."
1332 113 write (iout,*) "Error reading torsional energy parameters."
1334 114 write (iout,*) "Error reading double torsional energy parameters."
1337 & "Error reading cumulant (multibody energy) parameters."
1339 116 write (iout,*) "Error reading electrostatic energy parameters."
1341 1161 write (iout,*) "Error reading electrostatic energy parameters.Lip"
1343 117 write (iout,*) "Error reading side chain interaction parameters."
1345 118 write (iout,*) "Error reading SCp interaction parameters."
1347 119 write (iout,*) "Error reading SCCOR parameters"
1350 call MPI_Finalize(Ierror)
1357 subroutine getenv_loc(var, val)
1358 character(*) var, val
1361 character(2000) line
1364 open (196,file='env',status='old',readonly,shared)
1366 c write(*,*)'looking for ',var
1367 10 read(196,*,err=11,end=11)line
1368 iread=index(line,var)
1369 c write(*,*)iread,' ',var,' ',line
1370 if (iread.eq.0) go to 10
1371 c write(*,*)'---> ',line
1377 iread=iread+ilen(var)+1
1378 read (line(iread:),*,err=12,end=12) val
1379 c write(*,*)'OK: ',var,' = ',val
1385 #elif (defined CRAY)
1386 integer lennam,lenval,ierror
1388 c getenv using a POSIX call, useful on the T3D
1389 c Sept 1996, comment out error check on advice of H. Pritchard
1392 if(lennam.le.0) stop '--error calling getenv--'
1393 call pxfgetenv(var,lennam,val,lenval,ierror)
1394 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1396 call getenv(var,val)