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",ithep_pdb
394 read (ithep_pdb,*,err=111,end=111)
395 & a0thet(i),(athet(j,i,1,1),j=1,2),
396 & (bthet(j,i,1,1),j=1,2)
397 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
398 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
399 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
403 athet(1,i,1,-1)=athet(1,i,1,1)
404 athet(2,i,1,-1)=athet(2,i,1,1)
405 bthet(1,i,1,-1)=-bthet(1,i,1,1)
406 bthet(2,i,1,-1)=-bthet(2,i,1,1)
407 athet(1,i,-1,1)=-athet(1,i,1,1)
408 athet(2,i,-1,1)=-athet(2,i,1,1)
409 bthet(1,i,-1,1)=bthet(1,i,1,1)
410 bthet(2,i,-1,1)=bthet(2,i,1,1)
414 athet(1,i,-1,-1)=athet(1,-i,1,1)
415 athet(2,i,-1,-1)=-athet(2,-i,1,1)
416 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
417 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
418 athet(1,i,-1,1)=athet(1,-i,1,1)
419 athet(2,i,-1,1)=-athet(2,-i,1,1)
420 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
421 bthet(2,i,-1,1)=bthet(2,-i,1,1)
422 athet(1,i,1,-1)=-athet(1,-i,1,1)
423 athet(2,i,1,-1)=athet(2,-i,1,1)
424 bthet(1,i,1,-1)=bthet(1,-i,1,1)
425 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
430 polthet(j,i)=polthet(j,-i)
433 gthet(j,i)=gthet(j,-i)
436 write (2,*) "End reading THETA_PDB"
442 C Read the parameters of the probability distribution/energy expression
443 C of the side chains.
446 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
450 dsc_inv(i)=1.0D0/dsc(i)
461 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
462 & ((blower(k,l,1),l=1,k),k=1,3)
463 censc(1,1,-i)=censc(1,1,i)
464 censc(2,1,-i)=censc(2,1,i)
465 censc(3,1,-i)=-censc(3,1,i)
467 read (irotam,*,end=112,err=112) bsc(j,i)
468 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
469 & ((blower(k,l,j),l=1,k),k=1,3)
470 censc(1,j,-i)=censc(1,j,i)
471 censc(2,j,-i)=censc(2,j,i)
472 censc(3,j,-i)=-censc(3,j,i)
473 C BSC is amplitude of Gaussian
480 akl=akl+blower(k,m,j)*blower(l,m,j)
484 if (((k.eq.3).and.(l.ne.3))
485 & .or.((l.eq.3).and.(k.ne.3))) then
486 gaussc(k,l,j,-i)=-akl
487 gaussc(l,k,j,-i)=-akl
499 write (iout,'(/a)') 'Parameters of side-chain local geometry'
504 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
505 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
506 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
507 & 'log h',(bsc(j,i),j=1,nlobi)
508 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
509 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
511 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
512 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
515 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
516 write (iout,'(a,f10.4,4(16x,f10.4))')
517 & 'Center ',(bsc(j,i),j=1,nlobi)
518 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
527 C Read scrot parameters for potentials determined from all-atom AM1 calculations
528 C added by Urszula Kozlowska 07/11/2007
531 read (irotam,*,end=112,err=112)
533 read (irotam,*,end=112,err=112)
536 read(irotam,*,end=112,err=112) sc_parmin(j,i)
541 C Read the parameters of the probability distribution/energy expression
542 C of the side chains.
544 write (2,*) "Start reading ROTAM_PDB"
546 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
550 dsc_inv(i)=1.0D0/dsc(i)
561 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
562 & ((blower(k,l,1),l=1,k),k=1,3)
564 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
565 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
566 & ((blower(k,l,j),l=1,k),k=1,3)
573 akl=akl+blower(k,m,j)*blower(l,m,j)
583 write (2,*) "End reading ROTAM_PDB"
589 C Read torsional parameters in old format
591 read (itorp,*,end=113,err=113) ntortyp,nterm_old
592 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
593 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
598 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
604 write (iout,'(/a/)') 'Torsional constants:'
607 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
608 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
614 C Read torsional parameters
616 read (itorp,*,end=113,err=113) ntortyp
617 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
620 itortyp(i)=-itortyp(-i)
622 write (iout,*) 'ntortyp',ntortyp
624 do j=-ntortyp+1,ntortyp-1
625 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
627 nterm(-i,-j,iblock)=nterm(i,j,iblock)
628 nlor(-i,-j,iblock)=nlor(i,j,iblock)
631 do k=1,nterm(i,j,iblock)
632 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
634 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
635 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
636 v0ij=v0ij+si*v1(k,i,j,iblock)
638 c write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
639 c write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
640 c &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
642 do k=1,nlor(i,j,iblock)
643 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
644 & vlor2(k,i,j),vlor3(k,i,j)
645 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
648 v0(-i,-j,iblock)=v0ij
654 write (iout,'(/a/)') 'Torsional constants:'
657 write (iout,*) 'ityp',i,' jtyp',j
658 write (iout,*) 'Fourier constants'
659 do k=1,nterm(i,j,iblock)
660 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
663 write (iout,*) 'Lorenz constants'
664 do k=1,nlor(i,j,iblock)
665 write (iout,'(3(1pe15.5))')
666 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
673 C 6/23/01 Read parameters for double torsionals
677 do j=-ntortyp+1,ntortyp-1
678 do k=-ntortyp+1,ntortyp-1
679 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
680 c write (iout,*) "OK onelett",
683 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
684 & .or. t3.ne.toronelet(k)) then
685 write (iout,*) "Error in double torsional parameter file",
688 call MPI_Finalize(Ierror)
690 stop "Error in double torsional parameter file"
692 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
693 & ntermd_2(i,j,k,iblock)
694 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
695 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
696 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
697 & ntermd_1(i,j,k,iblock))
698 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
699 & ntermd_1(i,j,k,iblock))
700 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
701 & ntermd_1(i,j,k,iblock))
702 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
703 & ntermd_1(i,j,k,iblock))
704 C Martix of D parameters for one dimesional foureir series
705 do l=1,ntermd_1(i,j,k,iblock)
706 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
707 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
708 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
709 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
710 c write(iout,*) "whcodze" ,
711 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
713 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
714 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
715 & v2s(m,l,i,j,k,iblock),
716 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
717 C Martix of D parameters for two dimesional fourier series
718 do l=1,ntermd_2(i,j,k,iblock)
720 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
721 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
722 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
723 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
732 write (iout,*) 'Constants for double torsionals'
735 do j=-ntortyp+1,ntortyp-1
736 do k=-ntortyp+1,ntortyp-1
737 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
738 & ' nsingle',ntermd_1(i,j,k,iblock),
739 & ' ndouble',ntermd_2(i,j,k,iblock)
741 write (iout,*) 'Single angles:'
742 do l=1,ntermd_1(i,j,k,iblock)
743 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
744 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
745 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
746 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
749 write (iout,*) 'Pairs of angles:'
750 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
751 do l=1,ntermd_2(i,j,k,iblock)
752 write (iout,'(i5,20f10.5)')
753 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
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,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
760 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
769 C Read of Side-chain backbone correlation parameters
770 C Modified 11 May 2012 by Adasko
773 read (isccor,*,end=119,err=119) nsccortyp
775 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
777 isccortyp(i)=-isccortyp(-i)
779 iscprol=isccortyp(20)
780 c write (iout,*) 'ntortyp',ntortyp
782 cc maxinter is maximum interaction sites
786 read (isccor,*,end=119,err=119)
787 &nterm_sccor(i,j),nlor_sccor(i,j)
793 nterm_sccor(-i,j)=nterm_sccor(i,j)
794 nterm_sccor(-i,-j)=nterm_sccor(i,j)
795 nterm_sccor(i,-j)=nterm_sccor(i,j)
796 do k=1,nterm_sccor(i,j)
797 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
799 if (j.eq.iscprol) then
800 if (i.eq.isccortyp(10)) then
801 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
802 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
804 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
805 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
806 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
807 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
808 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
809 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
810 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
811 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
814 if (i.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 if (j.eq.isccortyp(10)) then
819 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
820 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)
824 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
825 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
826 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
827 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
831 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
832 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
833 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
834 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
837 do k=1,nlor_sccor(i,j)
838 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
839 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
840 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
841 &(1+vlor3sccor(k,i,j)**2)
843 v0sccor(l,i,j)=v0ijsccor
844 v0sccor(l,-i,j)=v0ijsccor1
845 v0sccor(l,i,-j)=v0ijsccor2
846 v0sccor(l,-i,-j)=v0ijsccor3
852 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
853 c write (iout,*) 'ntortyp',ntortyp
855 cc maxinter is maximum interaction sites
859 read (isccor,*,end=119,err=119)
860 & nterm_sccor(i,j),nlor_sccor(i,j)
864 do k=1,nterm_sccor(i,j)
865 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
867 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
870 do k=1,nlor_sccor(i,j)
871 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
872 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
873 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
874 &(1+vlor3sccor(k,i,j)**2)
876 v0sccor(i,j,iblock)=v0ijsccor
884 write (iout,'(/a/)') 'Torsional constants:'
887 write (iout,*) 'ityp',i,' jtyp',j
888 write (iout,*) 'Fourier constants'
889 do k=1,nterm_sccor(i,j)
890 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
892 write (iout,*) 'Lorenz constants'
893 do k=1,nlor_sccor(i,j)
894 write (iout,'(3(1pe15.5))')
895 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
902 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
903 C interaction energy of the Gly, Ala, and Pro prototypes.
907 write (iout,*) "Coefficients of the cumulants"
909 read (ifourier,*) nloctyp
912 read (ifourier,*,end=115,err=115)
913 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
915 read (ifourier,*,end=115,err=115) (bnew1(ii,1,i),ii=1,3)
916 read (ifourier,*,end=115,err=115) (bnew2(ii,1,i),ii=1,3)
917 read (ifourier,*,end=115,err=115) (bnew1(ii,2,i),ii=1,1)
918 read (ifourier,*,end=115,err=115) (bnew2(ii,2,i),ii=1,1)
919 read (ifourier,*,end=115,err=115) (eenew(ii,i),ii=1,1)
922 write (iout,*) 'Type',i
923 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
931 c B1tilde(1,i) = b(3)
932 c B1tilde(2,i) =-b(5)
933 c B1tilde(1,-i) =-b(3)
934 c B1tilde(2,-i) =b(5)
965 c Ctilde(1,1,i)=0.0d0
966 c Ctilde(1,2,i)=0.0d0
967 c Ctilde(2,1,i)=0.0d0
968 c Ctilde(2,2,i)=0.0d0
990 c Dtilde(1,1,i)=0.0d0
991 c Dtilde(1,2,i)=0.0d0
992 c Dtilde(2,1,i)=0.0d0
993 c Dtilde(2,2,i)=0.0d0
994 EE(1,1,i)= b(10)+b(11)
995 EE(2,2,i)=-b(10)+b(11)
996 EE(2,1,i)= b(12)-b(13)
997 EE(1,2,i)= b(12)+b(13)
998 EE(1,1,-i)= b(10)+b(11)
999 EE(2,2,-i)=-b(10)+b(11)
1000 EE(2,1,-i)=-b(12)+b(13)
1001 EE(1,2,-i)=-b(12)-b(13)
1007 c ee(2,1,i)=ee(1,2,i)
1011 write (iout,*) 'Type',i
1013 write(iout,*) B1(1,i),B1(2,i)
1015 write(iout,*) B2(1,i),B2(2,i)
1018 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1022 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1026 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1032 C Read electrostatic-interaction parameters
1036 write (iout,'(/a)') 'Electrostatic interaction constants:'
1037 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1038 & 'IT','JT','APP','BPP','AEL6','AEL3'
1040 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1041 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1042 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1043 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1048 app (i,j)=epp(i,j)*rri*rri
1049 bpp (i,j)=-2.0D0*epp(i,j)*rri
1050 ael6(i,j)=elpp6(i,j)*4.2D0**6
1051 ael3(i,j)=elpp3(i,j)*4.2D0**3
1053 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1054 & ael6(i,j),ael3(i,j)
1059 C Read side-chain interaction parameters.
1061 read (isidep,*,end=117,err=117) ipot,expon
1062 if (ipot.lt.1 .or. ipot.gt.5) then
1063 write (iout,'(2a)') 'Error while reading SC interaction',
1064 & 'potential file - unknown potential type.'
1066 call MPI_Finalize(Ierror)
1072 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1073 & ', exponents are ',expon,2*expon
1074 goto (10,20,30,30,40) ipot
1075 C----------------------- LJ potential ---------------------------------
1076 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1077 & (sigma0(i),i=1,ntyp)
1079 write (iout,'(/a/)') 'Parameters of the LJ 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,a)') 'residue','sigma'
1084 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1087 C----------------------- LJK potential --------------------------------
1088 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1089 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1091 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1092 write (iout,'(a/)') 'The epsilon array:'
1093 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1094 write (iout,'(/a)') 'One-body parameters:'
1095 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1096 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1100 C---------------------- GB or BP potential -----------------------------
1102 read (isidep,*,end=116,err=116)(eps(i,j),j=i,ntyp)
1104 read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
1105 read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
1106 read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
1107 read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
1108 C For the GB potential convert sigma'**2 into chi'
1111 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1115 write (iout,'(/a/)') 'Parameters of the BP potential:'
1116 write (iout,'(a/)') 'The epsilon array:'
1117 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1118 write (iout,'(/a)') 'One-body parameters:'
1119 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1121 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1122 & chip(i),alp(i),i=1,ntyp)
1125 C--------------------- GBV potential -----------------------------------
1126 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1127 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1128 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1130 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1131 write (iout,'(a/)') 'The epsilon array:'
1132 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1133 write (iout,'(/a)') 'One-body parameters:'
1134 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1135 & 's||/s_|_^2',' chip ',' alph '
1136 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1137 & sigii(i),chip(i),alp(i),i=1,ntyp)
1141 C-----------------------------------------------------------------------
1142 C Calculate the "working" parameters of SC interactions.
1150 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1151 sigma(j,i)=sigma(i,j)
1152 rs0(i,j)=dwa16*sigma(i,j)
1156 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1157 & 'Working parameters of the SC interactions:',
1158 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1163 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1172 sigeps=dsign(1.0D0,epsij)
1174 aa(i,j)=epsij*rrij*rrij
1175 bb(i,j)=-sigeps*epsij*rrij
1179 sigt1sq=sigma0(i)**2
1180 sigt2sq=sigma0(j)**2
1183 ratsig1=sigt2sq/sigt1sq
1184 ratsig2=1.0D0/ratsig1
1185 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1186 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1187 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1191 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1192 sigmaii(i,j)=rsum_max
1193 sigmaii(j,i)=rsum_max
1195 c sigmaii(i,j)=r0(i,j)
1196 c sigmaii(j,i)=r0(i,j)
1198 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1199 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1200 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1201 augm(i,j)=epsij*r_augm**(2*expon)
1202 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1209 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1210 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1211 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1217 C Define the SC-p interaction constants (hard-coded; old style)
1220 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1222 c aad(i,1)=0.3D0*4.0D0**12
1223 C Following line for constants currently implemented
1224 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1225 aad(i,1)=1.5D0*4.0D0**12
1226 c aad(i,1)=0.17D0*5.6D0**12
1228 C "Soft" SC-p repulsion
1230 C Following line for constants currently implemented
1231 c aad(i,1)=0.3D0*4.0D0**6
1232 C "Hard" SC-p repulsion
1233 bad(i,1)=3.0D0*4.0D0**6
1234 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1243 C 8/9/01 Read the SC-p interaction constants from file
1246 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1249 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1250 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1251 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1252 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1256 write (iout,*) "Parameters of SC-p interactions:"
1258 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1259 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1265 C Define the constants of the disulfide bridge
1269 c Old arbitrary potential - commented out.
1274 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1275 c energy surface of diethyl disulfide.
1276 c A. Liwo and U. Kozlowska, 11/24/03
1293 write (iout,'(/a)') "Disulfide bridge parameters:"
1294 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1295 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1296 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1297 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1301 111 write (iout,*) "Error reading bending energy parameters."
1303 112 write (iout,*) "Error reading rotamer energy parameters."
1305 113 write (iout,*) "Error reading torsional energy parameters."
1307 114 write (iout,*) "Error reading double torsional energy parameters."
1310 & "Error reading cumulant (multibody energy) parameters."
1312 116 write (iout,*) "Error reading electrostatic energy parameters."
1314 117 write (iout,*) "Error reading side chain interaction parameters."
1316 118 write (iout,*) "Error reading SCp interaction parameters."
1318 119 write (iout,*) "Error reading SCCOR parameters"
1321 call MPI_Finalize(Ierror)
1328 subroutine getenv_loc(var, val)
1329 character(*) var, val
1332 character(2000) line
1335 open (196,file='env',status='old',readonly,shared)
1337 c write(*,*)'looking for ',var
1338 10 read(196,*,err=11,end=11)line
1339 iread=index(line,var)
1340 c write(*,*)iread,' ',var,' ',line
1341 if (iread.eq.0) go to 10
1342 c write(*,*)'---> ',line
1348 iread=iread+ilen(var)+1
1349 read (line(iread:),*,err=12,end=12) val
1350 c write(*,*)'OK: ',var,' = ',val
1356 #elif (defined CRAY)
1357 integer lennam,lenval,ierror
1359 c getenv using a POSIX call, useful on the T3D
1360 c Sept 1996, comment out error check on advice of H. Pritchard
1363 if(lennam.le.0) stop '--error calling getenv--'
1364 call pxfgetenv(var,lennam,val,lenval,ierror)
1365 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1367 call getenv(var,val)