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,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,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)
102 C Read the parameters of the probability distribution/energy expression
103 C of the virtual-bond valence angles theta
106 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),
107 & (bthet(j,i,1,1),j=1,2)
108 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
109 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
110 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
114 athet(1,i,1,-1)=athet(1,i,1,1)
115 athet(2,i,1,-1)=athet(2,i,1,1)
116 bthet(1,i,1,-1)=-bthet(1,i,1,1)
117 bthet(2,i,1,-1)=-bthet(2,i,1,1)
118 athet(1,i,-1,1)=-athet(1,i,1,1)
119 athet(2,i,-1,1)=-athet(2,i,1,1)
120 bthet(1,i,-1,1)=bthet(1,i,1,1)
121 bthet(2,i,-1,1)=bthet(2,i,1,1)
125 athet(1,i,-1,-1)=athet(1,-i,1,1)
126 athet(2,i,-1,-1)=-athet(2,-i,1,1)
127 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
128 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
129 athet(1,i,-1,1)=athet(1,-i,1,1)
130 athet(2,i,-1,1)=-athet(2,-i,1,1)
131 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
132 bthet(2,i,-1,1)=bthet(2,-i,1,1)
133 athet(1,i,1,-1)=-athet(1,-i,1,1)
134 athet(2,i,1,-1)=athet(2,-i,1,1)
135 bthet(1,i,1,-1)=bthet(1,-i,1,1)
136 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
141 polthet(j,i)=polthet(j,-i)
144 gthet(j,i)=gthet(j,-i)
152 & 'Parameters of the virtual-bond valence angles:'
153 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
154 & ' ATHETA0 ',' A1 ',' A2 ',
157 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
158 & a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
160 write (iout,'(/a/9x,5a/79(1h-))')
161 & 'Parameters of the expression for sigma(theta_c):',
162 & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
163 & ' ALPH3 ',' SIGMA0C '
165 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
166 & (polthet(j,i),j=0,3),sigc0(i)
168 write (iout,'(/a/9x,5a/79(1h-))')
169 & 'Parameters of the second gaussian:',
170 & ' THETA0 ',' SIGMA0 ',' G1 ',
173 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
174 & sig0(i),(gthet(j,i),j=1,3)
178 & 'Parameters of the virtual-bond valence angles:'
179 write (iout,'(/a/9x,5a/79(1h-))')
180 & 'Coefficients of expansion',
181 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
182 & ' b1*10^1 ',' b2*10^1 '
184 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
185 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
186 & (10*bthet(j,i,1,1),j=1,2)
188 write (iout,'(/a/9x,5a/79(1h-))')
189 & 'Parameters of the expression for sigma(theta_c):',
190 & ' alpha0 ',' alph1 ',' alph2 ',
191 & ' alhp3 ',' sigma0c '
193 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
194 & (polthet(j,i),j=0,3),sigc0(i)
196 write (iout,'(/a/9x,5a/79(1h-))')
197 & 'Parameters of the second gaussian:',
198 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
201 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
202 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
208 C Read the parameters of Utheta determined from ab initio surfaces
209 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
211 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,
212 & ntheterm3,nsingle,ndouble
213 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
214 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
216 ithetyp(i)=-ithetyp(-i)
219 do i=-maxthetyp,maxthetyp
220 do j=-maxthetyp,maxthetyp
221 do k=-maxthetyp,maxthetyp
222 aa0thet(i,j,k,iblock)=0.0d0
224 aathet(l,i,j,k,iblock)=0.0d0
228 bbthet(m,l,i,j,k,iblock)=0.0d0
229 ccthet(m,l,i,j,k,iblock)=0.0d0
230 ddthet(m,l,i,j,k,iblock)=0.0d0
231 eethet(m,l,i,j,k,iblock)=0.0d0
237 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
238 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
246 c VAR:iblock means terminally blocking group 1=non-proline 2=proline
248 c VAR:ntethtyp is type of theta potentials type currently 0=glycine
249 c VAR:1=non-glicyne non-proline 2=proline
250 c VAR:negative values for D-aminoacid
252 do j=-nthetyp,nthetyp
253 do k=-nthetyp,nthetyp
254 read (ithep,'(6a)',end=111,err=111) res1
255 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
256 c VAR: aa0thet is variable describing the average value of Foureir
257 c VAR: expansion series
258 c VAR: aathet is foureir expansion in theta/2 angle for full formula
259 c VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
260 Condens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
261 read (ithep,*,end=111,err=111)
262 &(aathet(l,i,j,k,iblock),l=1,ntheterm)
263 read (ithep,*,end=111,err=111)
264 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
265 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
266 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
267 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),
269 read (ithep,*,end=111,err=111)
270 & (((ffthet(llll,lll,ll,i,j,k,iblock),
271 & ffthet(lll,llll,ll,i,j,k,iblock),
272 & ggthet(llll,lll,ll,i,j,k,iblock),
273 & ggthet(lll,llll,ll,i,j,k,iblock),
274 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
279 C For dummy ends assign glycine-type coefficients of theta-only terms; the
280 C coefficients of theta-and-gamma-dependent terms are zero.
281 C IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
282 C RECOMENTDED AFTER VERSION 3.3)
286 c aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
287 c aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
289 c aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
290 c aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
293 c aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
295 c aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
298 C AND COMMENT THE LOOPS BELOW
302 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
303 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
305 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
306 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
309 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
311 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
315 C Substitution for D aminoacids from symmetry.
318 do j=-nthetyp,nthetyp
319 do k=-nthetyp,nthetyp
320 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
322 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
326 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
327 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
328 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
329 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
335 ffthet(llll,lll,ll,i,j,k,iblock)=
336 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
337 ffthet(lll,llll,ll,i,j,k,iblock)=
338 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
339 ggthet(llll,lll,ll,i,j,k,iblock)=
340 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
341 ggthet(lll,llll,ll,i,j,k,iblock)=
342 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
351 C Control printout of the coefficients of virtual-bond-angle potentials
354 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
358 write (iout,'(//4a)')
359 & 'Type ',onelett(i),onelett(j),onelett(k)
360 write (iout,'(//a,10x,a)') " l","a[l]"
361 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
362 write (iout,'(i2,1pe15.5)')
363 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
365 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))')
366 & "b",l,"c",l,"d",l,"e",l
368 write (iout,'(i2,4(1pe15.5))') m,
369 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
370 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
374 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
375 & "f+",l,"f-",l,"g+",l,"g-",l
378 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
379 & ffthet(n,m,l,i,j,k,iblock),
380 & ffthet(m,n,l,i,j,k,iblock),
381 & ggthet(n,m,l,i,j,k,iblock),
382 & ggthet(m,n,l,i,j,k,iblock)
391 write (2,*) "Start reading THETA_PDB"
393 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),
394 & (bthet(j,i,1,1),j=1,2)
395 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
396 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
397 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
401 athet(1,i,1,-1)=athet(1,i,1,1)
402 athet(2,i,1,-1)=athet(2,i,1,1)
403 bthet(1,i,1,-1)=-bthet(1,i,1,1)
404 bthet(2,i,1,-1)=-bthet(2,i,1,1)
405 athet(1,i,-1,1)=-athet(1,i,1,1)
406 athet(2,i,-1,1)=-athet(2,i,1,1)
407 bthet(1,i,-1,1)=bthet(1,i,1,1)
408 bthet(2,i,-1,1)=bthet(2,i,1,1)
412 athet(1,i,-1,-1)=athet(1,-i,1,1)
413 athet(2,i,-1,-1)=-athet(2,-i,1,1)
414 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
415 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
416 athet(1,i,-1,1)=athet(1,-i,1,1)
417 athet(2,i,-1,1)=-athet(2,-i,1,1)
418 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
419 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)
428 polthet(j,i)=polthet(j,-i)
431 gthet(j,i)=gthet(j,-i)
434 write (2,*) "End reading THETA_PDB"
440 C Read the parameters of the probability distribution/energy expression
441 C of the side chains.
444 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
448 dsc_inv(i)=1.0D0/dsc(i)
459 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
460 & ((blower(k,l,1),l=1,k),k=1,3)
461 censc(1,1,-i)=censc(1,1,i)
462 censc(2,1,-i)=censc(2,1,i)
463 censc(3,1,-i)=-censc(3,1,i)
465 read (irotam,*,end=112,err=112) bsc(j,i)
466 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
467 & ((blower(k,l,j),l=1,k),k=1,3)
468 censc(1,j,-i)=censc(1,j,i)
469 censc(2,j,-i)=censc(2,j,i)
470 censc(3,j,-i)=-censc(3,j,i)
471 C BSC is amplitude of Gaussian
478 akl=akl+blower(k,m,j)*blower(l,m,j)
482 if (((k.eq.3).and.(l.ne.3))
483 & .or.((l.eq.3).and.(k.ne.3))) then
484 gaussc(k,l,j,-i)=-akl
485 gaussc(l,k,j,-i)=-akl
497 write (iout,'(/a)') 'Parameters of side-chain local geometry'
502 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
503 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
504 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
505 & 'log h',(bsc(j,i),j=1,nlobi)
506 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
507 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
509 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
510 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
513 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
514 write (iout,'(a,f10.4,4(16x,f10.4))')
515 & 'Center ',(bsc(j,i),j=1,nlobi)
516 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
525 C Read scrot parameters for potentials determined from all-atom AM1 calculations
526 C added by Urszula Kozlowska 07/11/2007
529 read (irotam,*,end=112,err=112)
531 read (irotam,*,end=112,err=112)
534 read(irotam,*,end=112,err=112) sc_parmin(j,i)
539 C Read the parameters of the probability distribution/energy expression
540 C of the side chains.
543 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
547 dsc_inv(i)=1.0D0/dsc(i)
558 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
559 & ((blower(k,l,1),l=1,k),k=1,3)
561 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
562 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
563 & ((blower(k,l,j),l=1,k),k=1,3)
570 akl=akl+blower(k,m,j)*blower(l,m,j)
585 C Read torsional parameters in old format
587 read (itorp,*,end=113,err=113) ntortyp,nterm_old
588 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
589 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
594 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
600 write (iout,'(/a/)') 'Torsional constants:'
603 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
604 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
610 C Read torsional parameters
612 read (itorp,*,end=113,err=113) ntortyp
613 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
616 itortyp(i)=-itortyp(-i)
618 write (iout,*) 'ntortyp',ntortyp
620 do j=-ntortyp+1,ntortyp-1
621 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
623 nterm(-i,-j,iblock)=nterm(i,j,iblock)
624 nlor(-i,-j,iblock)=nlor(i,j,iblock)
627 do k=1,nterm(i,j,iblock)
628 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
630 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
631 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
632 v0ij=v0ij+si*v1(k,i,j,iblock)
634 c write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
635 c write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
636 c &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
638 do k=1,nlor(i,j,iblock)
639 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
640 & vlor2(k,i,j),vlor3(k,i,j)
641 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
644 v0(-i,-j,iblock)=v0ij
650 write (iout,'(/a/)') 'Torsional constants:'
653 write (iout,*) 'ityp',i,' jtyp',j
654 write (iout,*) 'Fourier constants'
655 do k=1,nterm(i,j,iblock)
656 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
659 write (iout,*) 'Lorenz constants'
660 do k=1,nlor(i,j,iblock)
661 write (iout,'(3(1pe15.5))')
662 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
669 C 6/23/01 Read parameters for double torsionals
673 do j=-ntortyp+1,ntortyp-1
674 do k=-ntortyp+1,ntortyp-1
675 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
676 c write (iout,*) "OK onelett",
679 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
680 & .or. t3.ne.toronelet(k)) then
681 write (iout,*) "Error in double torsional parameter file",
684 call MPI_Finalize(Ierror)
686 stop "Error in double torsional parameter file"
688 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
689 & ntermd_2(i,j,k,iblock)
690 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
691 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
692 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
693 & ntermd_1(i,j,k,iblock))
694 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
695 & ntermd_1(i,j,k,iblock))
696 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
697 & ntermd_1(i,j,k,iblock))
698 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
699 & ntermd_1(i,j,k,iblock))
700 C Martix of D parameters for one dimesional foureir series
701 do l=1,ntermd_1(i,j,k,iblock)
702 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
703 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
704 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
705 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
706 c write(iout,*) "whcodze" ,
707 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
709 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
710 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
711 & v2s(m,l,i,j,k,iblock),
712 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
713 C Martix of D parameters for two dimesional fourier series
714 do l=1,ntermd_2(i,j,k,iblock)
716 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
717 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
718 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
719 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
728 write (iout,*) 'Constants for double torsionals'
731 do j=-ntortyp+1,ntortyp-1
732 do k=-ntortyp+1,ntortyp-1
733 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
734 & ' nsingle',ntermd_1(i,j,k,iblock),
735 & ' ndouble',ntermd_2(i,j,k,iblock)
737 write (iout,*) 'Single angles:'
738 do l=1,ntermd_1(i,j,k,iblock)
739 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
740 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
741 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
742 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
745 write (iout,*) 'Pairs of angles:'
746 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
747 do l=1,ntermd_2(i,j,k,iblock)
748 write (iout,'(i5,20f10.5)')
749 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
752 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
753 do l=1,ntermd_2(i,j,k,iblock)
754 write (iout,'(i5,20f10.5)')
755 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
756 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
765 C Read of Side-chain backbone correlation parameters
766 C Modified 11 May 2012 by Adasko
769 read (isccor,*,end=119,err=119) nsccortyp
771 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
773 isccortyp(i)=-isccortyp(-i)
775 iscprol=isccortyp(20)
776 c write (iout,*) 'ntortyp',ntortyp
778 cc maxinter is maximum interaction sites
782 read (isccor,*,end=119,err=119)
783 &nterm_sccor(i,j),nlor_sccor(i,j)
789 nterm_sccor(-i,j)=nterm_sccor(i,j)
790 nterm_sccor(-i,-j)=nterm_sccor(i,j)
791 nterm_sccor(i,-j)=nterm_sccor(i,j)
792 do k=1,nterm_sccor(i,j)
793 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
795 if (j.eq.iscprol) then
796 if (i.eq.isccortyp(10)) then
797 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
798 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
800 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
801 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
802 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
803 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
804 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
805 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
806 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
807 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
810 if (i.eq.isccortyp(10)) then
811 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
812 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
814 if (j.eq.isccortyp(10)) then
815 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
816 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
818 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
819 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
820 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
821 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
822 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
823 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
827 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
828 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
829 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
830 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
833 do k=1,nlor_sccor(i,j)
834 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
835 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
836 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
837 &(1+vlor3sccor(k,i,j)**2)
839 v0sccor(l,i,j)=v0ijsccor
840 v0sccor(l,-i,j)=v0ijsccor1
841 v0sccor(l,i,-j)=v0ijsccor2
842 v0sccor(l,-i,-j)=v0ijsccor3
848 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
849 c write (iout,*) 'ntortyp',ntortyp
851 cc maxinter is maximum interaction sites
855 read (isccor,*,end=119,err=119)
856 & nterm_sccor(i,j),nlor_sccor(i,j)
860 do k=1,nterm_sccor(i,j)
861 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
863 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
866 do k=1,nlor_sccor(i,j)
867 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
868 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
869 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
870 &(1+vlor3sccor(k,i,j)**2)
872 v0sccor(i,j,iblock)=v0ijsccor
880 write (iout,'(/a/)') 'Torsional constants:'
883 write (iout,*) 'ityp',i,' jtyp',j
884 write (iout,*) 'Fourier constants'
885 do k=1,nterm_sccor(i,j)
886 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
888 write (iout,*) 'Lorenz constants'
889 do k=1,nlor_sccor(i,j)
890 write (iout,'(3(1pe15.5))')
891 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
898 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
899 C interaction energy of the Gly, Ala, and Pro prototypes.
903 write (iout,*) "Coefficients of the cumulants"
905 read (ifourier,*) nloctyp
907 read (ifourier,*,end=115,err=115)
908 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
910 write (iout,*) 'Type',i
911 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
953 c Ctilde(1,1,i)=0.0d0
954 c Ctilde(1,2,i)=0.0d0
955 c Ctilde(2,1,i)=0.0d0
956 c Ctilde(2,2,i)=0.0d0
978 c Dtilde(1,1,i)=0.0d0
979 c Dtilde(1,2,i)=0.0d0
980 c Dtilde(2,1,i)=0.0d0
981 c Dtilde(2,2,i)=0.0d0
982 EE(1,1,i)= b(10)+b(11)
983 EE(2,2,i)=-b(10)+b(11)
984 EE(2,1,i)= b(12)-b(13)
985 EE(1,2,i)= b(12)+b(13)
986 EE(1,1,-i)= b(10)+b(11)
987 EE(2,2,-i)=-b(10)+b(11)
988 EE(2,1,-i)=-b(12)+b(13)
989 EE(1,2,-i)=-b(12)-b(13)
995 c ee(2,1,i)=ee(1,2,i)
999 write (iout,*) 'Type',i
1001 write(iout,*) B1(1,i),B1(2,i)
1003 write(iout,*) B2(1,i),B2(2,i)
1006 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1010 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1014 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1020 C Read electrostatic-interaction parameters
1024 write (iout,'(/a)') 'Electrostatic interaction constants:'
1025 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1026 & 'IT','JT','APP','BPP','AEL6','AEL3'
1028 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1029 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1030 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1031 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1036 app (i,j)=epp(i,j)*rri*rri
1037 bpp (i,j)=-2.0D0*epp(i,j)*rri
1038 ael6(i,j)=elpp6(i,j)*4.2D0**6
1039 ael3(i,j)=elpp3(i,j)*4.2D0**3
1041 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1042 & ael6(i,j),ael3(i,j)
1047 C Read side-chain interaction parameters.
1049 read (isidep,*,end=117,err=117) ipot,expon
1050 if (ipot.lt.1 .or. ipot.gt.5) then
1051 write (iout,'(2a)') 'Error while reading SC interaction',
1052 & 'potential file - unknown potential type.'
1054 call MPI_Finalize(Ierror)
1060 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1061 & ', exponents are ',expon,2*expon
1062 goto (10,20,30,30,40) ipot
1063 C----------------------- LJ potential ---------------------------------
1064 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1065 & (sigma0(i),i=1,ntyp)
1067 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1068 write (iout,'(a/)') 'The epsilon array:'
1069 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1070 write (iout,'(/a)') 'One-body parameters:'
1071 write (iout,'(a,4x,a)') 'residue','sigma'
1072 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1075 C----------------------- LJK potential --------------------------------
1076 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1077 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1079 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1080 write (iout,'(a/)') 'The epsilon array:'
1081 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1082 write (iout,'(/a)') 'One-body parameters:'
1083 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1084 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1088 C---------------------- GB or BP potential -----------------------------
1090 read (isidep,*,end=116,err=116)(eps(i,j),j=i,ntyp)
1092 read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
1093 read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
1094 read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
1095 read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
1096 C For the GB potential convert sigma'**2 into chi'
1099 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1103 write (iout,'(/a/)') 'Parameters of the BP potential:'
1104 write (iout,'(a/)') 'The epsilon array:'
1105 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1106 write (iout,'(/a)') 'One-body parameters:'
1107 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1109 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1110 & chip(i),alp(i),i=1,ntyp)
1113 C--------------------- GBV potential -----------------------------------
1114 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1115 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1116 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1118 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1119 write (iout,'(a/)') 'The epsilon array:'
1120 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1121 write (iout,'(/a)') 'One-body parameters:'
1122 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1123 & 's||/s_|_^2',' chip ',' alph '
1124 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1125 & sigii(i),chip(i),alp(i),i=1,ntyp)
1129 C-----------------------------------------------------------------------
1130 C Calculate the "working" parameters of SC interactions.
1138 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1139 sigma(j,i)=sigma(i,j)
1140 rs0(i,j)=dwa16*sigma(i,j)
1144 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1145 & 'Working parameters of the SC interactions:',
1146 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1151 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1160 sigeps=dsign(1.0D0,epsij)
1162 aa(i,j)=epsij*rrij*rrij
1163 bb(i,j)=-sigeps*epsij*rrij
1167 sigt1sq=sigma0(i)**2
1168 sigt2sq=sigma0(j)**2
1171 ratsig1=sigt2sq/sigt1sq
1172 ratsig2=1.0D0/ratsig1
1173 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1174 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1175 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1179 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1180 sigmaii(i,j)=rsum_max
1181 sigmaii(j,i)=rsum_max
1183 c sigmaii(i,j)=r0(i,j)
1184 c sigmaii(j,i)=r0(i,j)
1186 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1187 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1188 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1189 augm(i,j)=epsij*r_augm**(2*expon)
1190 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1197 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1198 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1199 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1205 C Define the SC-p interaction constants (hard-coded; old style)
1208 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1210 c aad(i,1)=0.3D0*4.0D0**12
1211 C Following line for constants currently implemented
1212 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1213 aad(i,1)=1.5D0*4.0D0**12
1214 c aad(i,1)=0.17D0*5.6D0**12
1216 C "Soft" SC-p repulsion
1218 C Following line for constants currently implemented
1219 c aad(i,1)=0.3D0*4.0D0**6
1220 C "Hard" SC-p repulsion
1221 bad(i,1)=3.0D0*4.0D0**6
1222 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1231 C 8/9/01 Read the SC-p interaction constants from file
1234 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1237 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1238 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1239 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1240 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1244 write (iout,*) "Parameters of SC-p interactions:"
1246 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1247 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1253 C Define the constants of the disulfide bridge
1257 c Old arbitrary potential - commented out.
1262 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1263 c energy surface of diethyl disulfide.
1264 c A. Liwo and U. Kozlowska, 11/24/03
1281 write (iout,'(/a)') "Disulfide bridge parameters:"
1282 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1283 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1284 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1285 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1289 111 write (iout,*) "Error reading bending energy parameters."
1291 112 write (iout,*) "Error reading rotamer energy parameters."
1293 113 write (iout,*) "Error reading torsional energy parameters."
1295 114 write (iout,*) "Error reading double torsional energy parameters."
1298 & "Error reading cumulant (multibody energy) parameters."
1300 116 write (iout,*) "Error reading electrostatic energy parameters."
1302 117 write (iout,*) "Error reading side chain interaction parameters."
1304 118 write (iout,*) "Error reading SCp interaction parameters."
1306 119 write (iout,*) "Error reading SCCOR parameters"
1309 call MPI_Finalize(Ierror)
1316 subroutine getenv_loc(var, val)
1317 character(*) var, val
1320 character(2000) line
1323 open (196,file='env',status='old',readonly,shared)
1325 c write(*,*)'looking for ',var
1326 10 read(196,*,err=11,end=11)line
1327 iread=index(line,var)
1328 c write(*,*)iread,' ',var,' ',line
1329 if (iread.eq.0) go to 10
1330 c write(*,*)'---> ',line
1336 iread=iread+ilen(var)+1
1337 read (line(iread:),*,err=12,end=12) val
1338 c write(*,*)'OK: ',var,' = ',val
1344 #elif (defined CRAY)
1345 integer lennam,lenval,ierror
1347 c getenv using a POSIX call, useful on the T3D
1348 c Sept 1996, comment out error check on advice of H. Pritchard
1351 if(lennam.le.0) stop '--error calling getenv--'
1352 call pxfgetenv(var,lennam,val,lenval,ierror)
1353 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1355 call getenv(var,val)