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
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)
1008 write(iout,*) "TU DOCHODZE"
1014 c ee(2,1,i)=ee(1,2,i)
1019 write (iout,*) 'Type',i
1021 write(iout,*) B1(1,i),B1(2,i)
1023 write(iout,*) B2(1,i),B2(2,i)
1026 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1030 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1034 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1041 C Read electrostatic-interaction parameters
1045 write (iout,'(/a)') 'Electrostatic interaction constants:'
1046 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1047 & 'IT','JT','APP','BPP','AEL6','AEL3'
1049 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1050 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1051 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1052 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1057 app (i,j)=epp(i,j)*rri*rri
1058 bpp (i,j)=-2.0D0*epp(i,j)*rri
1059 ael6(i,j)=elpp6(i,j)*4.2D0**6
1060 ael3(i,j)=elpp3(i,j)*4.2D0**3
1062 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1063 & ael6(i,j),ael3(i,j)
1068 C Read side-chain interaction parameters.
1070 read (isidep,*,end=117,err=117) ipot,expon
1071 if (ipot.lt.1 .or. ipot.gt.5) then
1072 write (iout,'(2a)') 'Error while reading SC interaction',
1073 & 'potential file - unknown potential type.'
1075 call MPI_Finalize(Ierror)
1081 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1082 & ', exponents are ',expon,2*expon
1083 goto (10,20,30,30,40) ipot
1084 C----------------------- LJ potential ---------------------------------
1085 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1086 & (sigma0(i),i=1,ntyp)
1088 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1089 write (iout,'(a/)') 'The epsilon array:'
1090 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1091 write (iout,'(/a)') 'One-body parameters:'
1092 write (iout,'(a,4x,a)') 'residue','sigma'
1093 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1096 C----------------------- LJK potential --------------------------------
1097 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1098 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1100 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1101 write (iout,'(a/)') 'The epsilon array:'
1102 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1103 write (iout,'(/a)') 'One-body parameters:'
1104 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1105 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1109 C---------------------- GB or BP potential -----------------------------
1111 read (isidep,*,end=116,err=116)(eps(i,j),j=i,ntyp)
1113 read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
1114 read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
1115 read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
1116 read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
1117 C now we start reading lipid
1119 read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
1120 C print *,"WARNING!!"
1122 C epslip(i,j)=epslip(i,j)+0.05d0
1125 write(iout,*) epslip(1,1),"OK?"
1126 C For the GB potential convert sigma'**2 into chi'
1129 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1133 write (iout,'(/a/)') 'Parameters of the BP potential:'
1134 write (iout,'(a/)') 'The epsilon array:'
1135 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1136 write (iout,'(/a)') 'One-body parameters:'
1137 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1139 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1140 & chip(i),alp(i),i=1,ntyp)
1143 C--------------------- GBV potential -----------------------------------
1144 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1145 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1146 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1148 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1149 write (iout,'(a/)') 'The epsilon array:'
1150 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1151 write (iout,'(/a)') 'One-body parameters:'
1152 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1153 & 's||/s_|_^2',' chip ',' alph '
1154 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1155 & sigii(i),chip(i),alp(i),i=1,ntyp)
1159 C-----------------------------------------------------------------------
1160 C Calculate the "working" parameters of SC interactions.
1164 epslip(i,j)=epslip(j,i)
1169 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1170 sigma(j,i)=sigma(i,j)
1171 rs0(i,j)=dwa16*sigma(i,j)
1175 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1176 & 'Working parameters of the SC interactions:',
1177 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1182 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1191 sigeps=dsign(1.0D0,epsij)
1193 aa_aq(i,j)=epsij*rrij*rrij
1194 bb_aq(i,j)=-sigeps*epsij*rrij
1195 aa_aq(j,i)=aa_aq(i,j)
1196 bb_aq(j,i)=bb_aq(i,j)
1197 epsijlip=epslip(i,j)
1198 sigeps=dsign(1.0D0,epsijlip)
1199 epsijlip=dabs(epsijlip)
1200 aa_lip(i,j)=epsijlip*rrij*rrij
1201 bb_lip(i,j)=-sigeps*epsijlip*rrij
1202 aa_lip(j,i)=aa_lip(i,j)
1203 bb_lip(j,i)=bb_lip(i,j)
1205 sigt1sq=sigma0(i)**2
1206 sigt2sq=sigma0(j)**2
1209 ratsig1=sigt2sq/sigt1sq
1210 ratsig2=1.0D0/ratsig1
1211 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1212 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1213 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1217 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1218 sigmaii(i,j)=rsum_max
1219 sigmaii(j,i)=rsum_max
1221 c sigmaii(i,j)=r0(i,j)
1222 c sigmaii(j,i)=r0(i,j)
1224 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1225 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1226 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1227 augm(i,j)=epsij*r_augm**(2*expon)
1228 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1235 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1236 & restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1237 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1243 C Define the SC-p interaction constants (hard-coded; old style)
1246 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1248 c aad(i,1)=0.3D0*4.0D0**12
1249 C Following line for constants currently implemented
1250 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1251 aad(i,1)=1.5D0*4.0D0**12
1252 c aad(i,1)=0.17D0*5.6D0**12
1254 C "Soft" SC-p repulsion
1256 C Following line for constants currently implemented
1257 c aad(i,1)=0.3D0*4.0D0**6
1258 C "Hard" SC-p repulsion
1259 bad(i,1)=3.0D0*4.0D0**6
1260 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1269 C 8/9/01 Read the SC-p interaction constants from file
1272 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1275 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1276 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1277 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1278 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1282 write (iout,*) "Parameters of SC-p interactions:"
1284 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1285 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1291 C Define the constants of the disulfide bridge
1295 c Old arbitrary potential - commented out.
1300 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1301 c energy surface of diethyl disulfide.
1302 c A. Liwo and U. Kozlowska, 11/24/03
1319 write (iout,'(/a)') "Disulfide bridge parameters:"
1320 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1321 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1322 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1323 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1327 111 write (iout,*) "Error reading bending energy parameters."
1329 112 write (iout,*) "Error reading rotamer energy parameters."
1331 113 write (iout,*) "Error reading torsional energy parameters."
1333 114 write (iout,*) "Error reading double torsional energy parameters."
1336 & "Error reading cumulant (multibody energy) parameters."
1338 116 write (iout,*) "Error reading electrostatic energy parameters."
1340 1161 write (iout,*) "Error reading electrostatic energy parameters.Lip"
1342 117 write (iout,*) "Error reading side chain interaction parameters."
1344 118 write (iout,*) "Error reading SCp interaction parameters."
1346 119 write (iout,*) "Error reading SCCOR parameters"
1349 call MPI_Finalize(Ierror)
1356 subroutine getenv_loc(var, val)
1357 character(*) var, val
1360 character(2000) line
1363 open (196,file='env',status='old',readonly,shared)
1365 c write(*,*)'looking for ',var
1366 10 read(196,*,err=11,end=11)line
1367 iread=index(line,var)
1368 c write(*,*)iread,' ',var,' ',line
1369 if (iread.eq.0) go to 10
1370 c write(*,*)'---> ',line
1376 iread=iread+ilen(var)+1
1377 read (line(iread:),*,err=12,end=12) val
1378 c write(*,*)'OK: ',var,' = ',val
1384 #elif (defined CRAY)
1385 integer lennam,lenval,ierror
1387 c getenv using a POSIX call, useful on the T3D
1388 c Sept 1996, comment out error check on advice of H. Pritchard
1391 if(lennam.le.0) stop '--error calling getenv--'
1392 call pxfgetenv(var,lennam,val,lenval,ierror)
1393 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1395 call getenv(var,val)