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
911 read (ifourier,*,end=115,err=115)
912 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
914 write (iout,*) 'Type',i
915 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
957 c Ctilde(1,1,i)=0.0d0
958 c Ctilde(1,2,i)=0.0d0
959 c Ctilde(2,1,i)=0.0d0
960 c Ctilde(2,2,i)=0.0d0
982 c Dtilde(1,1,i)=0.0d0
983 c Dtilde(1,2,i)=0.0d0
984 c Dtilde(2,1,i)=0.0d0
985 c Dtilde(2,2,i)=0.0d0
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)
990 EE(1,1,-i)= b(10)+b(11)
991 EE(2,2,-i)=-b(10)+b(11)
992 EE(2,1,-i)=-b(12)+b(13)
993 EE(1,2,-i)=-b(12)-b(13)
999 c ee(2,1,i)=ee(1,2,i)
1003 write (iout,*) 'Type',i
1005 write(iout,*) B1(1,i),B1(2,i)
1007 write(iout,*) B2(1,i),B2(2,i)
1010 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1014 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1018 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1024 C Read electrostatic-interaction parameters
1028 write (iout,'(/a)') 'Electrostatic interaction constants:'
1029 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1030 & 'IT','JT','APP','BPP','AEL6','AEL3'
1032 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1033 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1034 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1035 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1040 app (i,j)=epp(i,j)*rri*rri
1041 bpp (i,j)=-2.0D0*epp(i,j)*rri
1042 ael6(i,j)=elpp6(i,j)*4.2D0**6
1043 ael3(i,j)=elpp3(i,j)*4.2D0**3
1045 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1046 & ael6(i,j),ael3(i,j)
1051 C Read side-chain interaction parameters.
1053 read (isidep,*,end=117,err=117) ipot,expon
1054 if (ipot.lt.1 .or. ipot.gt.5) then
1055 write (iout,'(2a)') 'Error while reading SC interaction',
1056 & 'potential file - unknown potential type.'
1058 call MPI_Finalize(Ierror)
1064 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1065 & ', exponents are ',expon,2*expon
1066 goto (10,20,30,30,40) ipot
1067 C----------------------- LJ potential ---------------------------------
1068 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1069 & (sigma0(i),i=1,ntyp)
1071 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1072 write (iout,'(a/)') 'The epsilon array:'
1073 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1074 write (iout,'(/a)') 'One-body parameters:'
1075 write (iout,'(a,4x,a)') 'residue','sigma'
1076 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1079 C----------------------- LJK potential --------------------------------
1080 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1081 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1083 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1084 write (iout,'(a/)') 'The epsilon array:'
1085 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1086 write (iout,'(/a)') 'One-body parameters:'
1087 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1088 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1092 C---------------------- GB or BP potential -----------------------------
1094 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1096 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
1097 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
1098 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
1099 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
1100 C For the GB potential convert sigma'**2 into chi'
1103 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1107 write (iout,'(/a/)') 'Parameters of the BP potential:'
1108 write (iout,'(a/)') 'The epsilon array:'
1109 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1110 write (iout,'(/a)') 'One-body parameters:'
1111 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1113 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1114 & chip(i),alp(i),i=1,ntyp)
1117 C--------------------- GBV potential -----------------------------------
1118 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1119 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1120 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1122 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1123 write (iout,'(a/)') 'The epsilon array:'
1124 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1125 write (iout,'(/a)') 'One-body parameters:'
1126 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1127 & 's||/s_|_^2',' chip ',' alph '
1128 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1129 & sigii(i),chip(i),alp(i),i=1,ntyp)
1133 C-----------------------------------------------------------------------
1134 C Calculate the "working" parameters of SC interactions.
1142 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1143 sigma(j,i)=sigma(i,j)
1144 rs0(i,j)=dwa16*sigma(i,j)
1148 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1149 & 'Working parameters of the SC interactions:',
1150 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1155 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1164 sigeps=dsign(1.0D0,epsij)
1166 aa(i,j)=epsij*rrij*rrij
1167 bb(i,j)=-sigeps*epsij*rrij
1171 sigt1sq=sigma0(i)**2
1172 sigt2sq=sigma0(j)**2
1175 ratsig1=sigt2sq/sigt1sq
1176 ratsig2=1.0D0/ratsig1
1177 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1178 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1179 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1183 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1184 sigmaii(i,j)=rsum_max
1185 sigmaii(j,i)=rsum_max
1187 c sigmaii(i,j)=r0(i,j)
1188 c sigmaii(j,i)=r0(i,j)
1190 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1191 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1192 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1193 augm(i,j)=epsij*r_augm**(2*expon)
1194 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1201 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1202 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1203 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1209 C Define the SC-p interaction constants (hard-coded; old style)
1212 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1214 c aad(i,1)=0.3D0*4.0D0**12
1215 C Following line for constants currently implemented
1216 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1217 aad(i,1)=1.5D0*4.0D0**12
1218 c aad(i,1)=0.17D0*5.6D0**12
1220 C "Soft" SC-p repulsion
1222 C Following line for constants currently implemented
1223 c aad(i,1)=0.3D0*4.0D0**6
1224 C "Hard" SC-p repulsion
1225 bad(i,1)=3.0D0*4.0D0**6
1226 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1235 C 8/9/01 Read the SC-p interaction constants from file
1238 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1241 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1242 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1243 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1244 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1248 write (iout,*) "Parameters of SC-p interactions:"
1250 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1251 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1257 C Define the constants of the disulfide bridge
1261 c Old arbitrary potential - commented out.
1266 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1267 c energy surface of diethyl disulfide.
1268 c A. Liwo and U. Kozlowska, 11/24/03
1285 write (iout,'(/a)') "Disulfide bridge parameters:"
1286 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1287 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1288 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1289 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1293 111 write (iout,*) "Error reading bending energy parameters."
1295 112 write (iout,*) "Error reading rotamer energy parameters."
1297 113 write (iout,*) "Error reading torsional energy parameters."
1299 114 write (iout,*) "Error reading double torsional energy parameters."
1302 & "Error reading cumulant (multibody energy) parameters."
1304 116 write (iout,*) "Error reading electrostatic energy parameters."
1306 117 write (iout,*) "Error reading side chain interaction parameters."
1308 118 write (iout,*) "Error reading SCp interaction parameters."
1310 119 write (iout,*) "Error reading SCCOR parameters"
1313 call MPI_Finalize(Ierror)
1320 subroutine getenv_loc(var, val)
1321 character(*) var, val
1324 character(2000) line
1327 open (196,file='env',status='old',readonly,shared)
1329 c write(*,*)'looking for ',var
1330 10 read(196,*,err=11,end=11)line
1331 iread=index(line,var)
1332 c write(*,*)iread,' ',var,' ',line
1333 if (iread.eq.0) go to 10
1334 c write(*,*)'---> ',line
1340 iread=iread+ilen(var)+1
1341 read (line(iread:),*,err=12,end=12) val
1342 c write(*,*)'OK: ',var,' = ',val
1348 #elif (defined CRAY)
1349 integer lennam,lenval,ierror
1351 c getenv using a POSIX call, useful on the T3D
1352 c Sept 1996, comment out error check on advice of H. Pritchard
1355 if(lennam.le.0) stop '--error calling getenv--'
1356 call pxfgetenv(var,lennam,val,lenval,ierror)
1357 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1359 call getenv(var,val)