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
77 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
78 & j=1,nbondterm(i)),msc(i),isc(i),restok(i)
83 dsc_inv(i)=1.0D0/dsc(i)
88 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
89 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',
91 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
93 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
94 & vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
96 write (iout,'(13x,3f10.5)')
97 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
103 C Read the parameters of the probability distribution/energy expression
104 C of the virtual-bond valence angles theta
107 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),
108 & (bthet(j,i,1,1),j=1,2)
109 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
110 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
111 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
115 athet(1,i,1,-1)=athet(1,i,1,1)
116 athet(2,i,1,-1)=athet(2,i,1,1)
117 bthet(1,i,1,-1)=-bthet(1,i,1,1)
118 bthet(2,i,1,-1)=-bthet(2,i,1,1)
119 athet(1,i,-1,1)=-athet(1,i,1,1)
120 athet(2,i,-1,1)=-athet(2,i,1,1)
121 bthet(1,i,-1,1)=bthet(1,i,1,1)
122 bthet(2,i,-1,1)=bthet(2,i,1,1)
126 athet(1,i,-1,-1)=athet(1,-i,1,1)
127 athet(2,i,-1,-1)=-athet(2,-i,1,1)
128 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
129 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
130 athet(1,i,-1,1)=athet(1,-i,1,1)
131 athet(2,i,-1,1)=-athet(2,-i,1,1)
132 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
133 bthet(2,i,-1,1)=bthet(2,-i,1,1)
134 athet(1,i,1,-1)=-athet(1,-i,1,1)
135 athet(2,i,1,-1)=athet(2,-i,1,1)
136 bthet(1,i,1,-1)=bthet(1,-i,1,1)
137 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
142 polthet(j,i)=polthet(j,-i)
145 gthet(j,i)=gthet(j,-i)
153 & 'Parameters of the virtual-bond valence angles:'
154 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
155 & ' ATHETA0 ',' A1 ',' A2 ',
158 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
159 & a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
161 write (iout,'(/a/9x,5a/79(1h-))')
162 & 'Parameters of the expression for sigma(theta_c):',
163 & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
164 & ' ALPH3 ',' SIGMA0C '
166 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
167 & (polthet(j,i),j=0,3),sigc0(i)
169 write (iout,'(/a/9x,5a/79(1h-))')
170 & 'Parameters of the second gaussian:',
171 & ' THETA0 ',' SIGMA0 ',' G1 ',
174 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
175 & sig0(i),(gthet(j,i),j=1,3)
179 & 'Parameters of the virtual-bond valence angles:'
180 write (iout,'(/a/9x,5a/79(1h-))')
181 & 'Coefficients of expansion',
182 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
183 & ' b1*10^1 ',' b2*10^1 '
185 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
186 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
187 & (10*bthet(j,i,1,1),j=1,2)
189 write (iout,'(/a/9x,5a/79(1h-))')
190 & 'Parameters of the expression for sigma(theta_c):',
191 & ' alpha0 ',' alph1 ',' alph2 ',
192 & ' alhp3 ',' sigma0c '
194 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
195 & (polthet(j,i),j=0,3),sigc0(i)
197 write (iout,'(/a/9x,5a/79(1h-))')
198 & 'Parameters of the second gaussian:',
199 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
202 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
203 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
209 C Read the parameters of Utheta determined from ab initio surfaces
210 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
212 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,
213 & ntheterm3,nsingle,ndouble
214 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
215 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
217 ithetyp(i)=-ithetyp(-i)
220 do i=-maxthetyp,maxthetyp
221 do j=-maxthetyp,maxthetyp
222 do k=-maxthetyp,maxthetyp
223 aa0thet(i,j,k,iblock)=0.0d0
225 aathet(l,i,j,k,iblock)=0.0d0
229 bbthet(m,l,i,j,k,iblock)=0.0d0
230 ccthet(m,l,i,j,k,iblock)=0.0d0
231 ddthet(m,l,i,j,k,iblock)=0.0d0
232 eethet(m,l,i,j,k,iblock)=0.0d0
238 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
239 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
247 c VAR:iblock means terminally blocking group 1=non-proline 2=proline
249 c VAR:ntethtyp is type of theta potentials type currently 0=glycine
250 c VAR:1=non-glicyne non-proline 2=proline
251 c VAR:negative values for D-aminoacid
253 do j=-nthetyp,nthetyp
254 do k=-nthetyp,nthetyp
255 read (ithep,'(6a)',end=111,err=111) res1
256 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
257 c VAR: aa0thet is variable describing the average value of Foureir
258 c VAR: expansion series
259 c VAR: aathet is foureir expansion in theta/2 angle for full formula
260 c VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
261 Condens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
262 read (ithep,*,end=111,err=111)
263 &(aathet(l,i,j,k,iblock),l=1,ntheterm)
264 read (ithep,*,end=111,err=111)
265 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
266 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
267 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
268 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),
270 read (ithep,*,end=111,err=111)
271 & (((ffthet(llll,lll,ll,i,j,k,iblock),
272 & ffthet(lll,llll,ll,i,j,k,iblock),
273 & ggthet(llll,lll,ll,i,j,k,iblock),
274 & ggthet(lll,llll,ll,i,j,k,iblock),
275 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
280 C For dummy ends assign glycine-type coefficients of theta-only terms; the
281 C coefficients of theta-and-gamma-dependent terms are zero.
282 C IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
283 C RECOMENTDED AFTER VERSION 3.3)
287 c aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
288 c aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
290 c aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
291 c aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
294 c aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
296 c aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
299 C AND COMMENT THE LOOPS BELOW
303 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
304 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
306 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
307 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
310 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
312 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
316 C Substitution for D aminoacids from symmetry.
319 do j=-nthetyp,nthetyp
320 do k=-nthetyp,nthetyp
321 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
323 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
327 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
328 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
329 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
330 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
336 ffthet(llll,lll,ll,i,j,k,iblock)=
337 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
338 ffthet(lll,llll,ll,i,j,k,iblock)=
339 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
340 ggthet(llll,lll,ll,i,j,k,iblock)=
341 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
342 ggthet(lll,llll,ll,i,j,k,iblock)=
343 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
352 C Control printout of the coefficients of virtual-bond-angle potentials
355 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
359 write (iout,'(//4a)')
360 & 'Type ',onelett(i),onelett(j),onelett(k)
361 write (iout,'(//a,10x,a)') " l","a[l]"
362 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
363 write (iout,'(i2,1pe15.5)')
364 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
366 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))')
367 & "b",l,"c",l,"d",l,"e",l
369 write (iout,'(i2,4(1pe15.5))') m,
370 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
371 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
375 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
376 & "f+",l,"f-",l,"g+",l,"g-",l
379 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
380 & ffthet(n,m,l,i,j,k,iblock),
381 & ffthet(m,n,l,i,j,k,iblock),
382 & ggthet(n,m,l,i,j,k,iblock),
383 & ggthet(m,n,l,i,j,k,iblock)
392 write (2,*) "Start reading THETA_PDB",ithep_pdb
395 read (ithep_pdb,*,err=111,end=111)
396 & a0thet(i),(athet(j,i,1,1),j=1,2),
397 & (bthet(j,i,1,1),j=1,2)
398 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
399 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
400 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
404 athet(1,i,1,-1)=athet(1,i,1,1)
405 athet(2,i,1,-1)=athet(2,i,1,1)
406 bthet(1,i,1,-1)=-bthet(1,i,1,1)
407 bthet(2,i,1,-1)=-bthet(2,i,1,1)
408 athet(1,i,-1,1)=-athet(1,i,1,1)
409 athet(2,i,-1,1)=-athet(2,i,1,1)
410 bthet(1,i,-1,1)=bthet(1,i,1,1)
411 bthet(2,i,-1,1)=bthet(2,i,1,1)
415 athet(1,i,-1,-1)=athet(1,-i,1,1)
416 athet(2,i,-1,-1)=-athet(2,-i,1,1)
417 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
418 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
419 athet(1,i,-1,1)=athet(1,-i,1,1)
420 athet(2,i,-1,1)=-athet(2,-i,1,1)
421 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
422 bthet(2,i,-1,1)=bthet(2,-i,1,1)
423 athet(1,i,1,-1)=-athet(1,-i,1,1)
424 athet(2,i,1,-1)=athet(2,-i,1,1)
425 bthet(1,i,1,-1)=bthet(1,-i,1,1)
426 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
431 polthet(j,i)=polthet(j,-i)
434 gthet(j,i)=gthet(j,-i)
437 write (2,*) "End reading THETA_PDB"
443 C Read the parameters of the probability distribution/energy expression
444 C of the side chains.
447 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
451 dsc_inv(i)=1.0D0/dsc(i)
462 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
463 & ((blower(k,l,1),l=1,k),k=1,3)
464 censc(1,1,-i)=censc(1,1,i)
465 censc(2,1,-i)=censc(2,1,i)
466 censc(3,1,-i)=-censc(3,1,i)
468 read (irotam,*,end=112,err=112) bsc(j,i)
469 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
470 & ((blower(k,l,j),l=1,k),k=1,3)
471 censc(1,j,-i)=censc(1,j,i)
472 censc(2,j,-i)=censc(2,j,i)
473 censc(3,j,-i)=-censc(3,j,i)
474 C BSC is amplitude of Gaussian
481 akl=akl+blower(k,m,j)*blower(l,m,j)
485 if (((k.eq.3).and.(l.ne.3))
486 & .or.((l.eq.3).and.(k.ne.3))) then
487 gaussc(k,l,j,-i)=-akl
488 gaussc(l,k,j,-i)=-akl
500 write (iout,'(/a)') 'Parameters of side-chain local geometry'
505 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
506 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
507 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
508 & 'log h',(bsc(j,i),j=1,nlobi)
509 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
510 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
512 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
513 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
516 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
517 write (iout,'(a,f10.4,4(16x,f10.4))')
518 & 'Center ',(bsc(j,i),j=1,nlobi)
519 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
528 C Read scrot parameters for potentials determined from all-atom AM1 calculations
529 C added by Urszula Kozlowska 07/11/2007
532 read (irotam,*,end=112,err=112)
534 read (irotam,*,end=112,err=112)
537 read(irotam,*,end=112,err=112) sc_parmin(j,i)
542 C Read the parameters of the probability distribution/energy expression
543 C of the side chains.
545 write (2,*) "Start reading ROTAM_PDB"
547 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
551 dsc_inv(i)=1.0D0/dsc(i)
562 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
563 & ((blower(k,l,1),l=1,k),k=1,3)
565 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
566 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
567 & ((blower(k,l,j),l=1,k),k=1,3)
574 akl=akl+blower(k,m,j)*blower(l,m,j)
584 write (2,*) "End reading ROTAM_PDB"
590 C Read torsional parameters in old format
592 read (itorp,*,end=113,err=113) ntortyp,nterm_old
593 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
594 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
599 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
605 write (iout,'(/a/)') 'Torsional constants:'
608 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
609 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
615 C Read torsional parameters
617 read (itorp,*,end=113,err=113) ntortyp
618 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
621 itortyp(i)=-itortyp(-i)
623 write (iout,*) 'ntortyp',ntortyp
625 do j=-ntortyp+1,ntortyp-1
626 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
628 nterm(-i,-j,iblock)=nterm(i,j,iblock)
629 nlor(-i,-j,iblock)=nlor(i,j,iblock)
632 do k=1,nterm(i,j,iblock)
633 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
635 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
636 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
637 v0ij=v0ij+si*v1(k,i,j,iblock)
639 c write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
640 c write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
641 c &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
643 do k=1,nlor(i,j,iblock)
644 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
645 & vlor2(k,i,j),vlor3(k,i,j)
646 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
649 v0(-i,-j,iblock)=v0ij
655 write (iout,'(/a/)') 'Torsional constants:'
658 write (iout,*) 'ityp',i,' jtyp',j
659 write (iout,*) 'Fourier constants'
660 do k=1,nterm(i,j,iblock)
661 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
664 write (iout,*) 'Lorenz constants'
665 do k=1,nlor(i,j,iblock)
666 write (iout,'(3(1pe15.5))')
667 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
674 C 6/23/01 Read parameters for double torsionals
678 do j=-ntortyp+1,ntortyp-1
679 do k=-ntortyp+1,ntortyp-1
680 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
681 c write (iout,*) "OK onelett",
684 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
685 & .or. t3.ne.toronelet(k)) then
686 write (iout,*) "Error in double torsional parameter file",
689 call MPI_Finalize(Ierror)
691 stop "Error in double torsional parameter file"
693 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
694 & ntermd_2(i,j,k,iblock)
695 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
696 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
697 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
698 & ntermd_1(i,j,k,iblock))
699 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
700 & ntermd_1(i,j,k,iblock))
701 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
702 & ntermd_1(i,j,k,iblock))
703 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
704 & ntermd_1(i,j,k,iblock))
705 C Martix of D parameters for one dimesional foureir series
706 do l=1,ntermd_1(i,j,k,iblock)
707 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
708 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
709 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
710 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
711 c write(iout,*) "whcodze" ,
712 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
714 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
715 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
716 & v2s(m,l,i,j,k,iblock),
717 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
718 C Martix of D parameters for two dimesional fourier series
719 do l=1,ntermd_2(i,j,k,iblock)
721 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
722 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
723 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
724 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
733 write (iout,*) 'Constants for double torsionals'
736 do j=-ntortyp+1,ntortyp-1
737 do k=-ntortyp+1,ntortyp-1
738 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
739 & ' nsingle',ntermd_1(i,j,k,iblock),
740 & ' ndouble',ntermd_2(i,j,k,iblock)
742 write (iout,*) 'Single angles:'
743 do l=1,ntermd_1(i,j,k,iblock)
744 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
745 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
746 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
747 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
750 write (iout,*) 'Pairs of angles:'
751 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
752 do l=1,ntermd_2(i,j,k,iblock)
753 write (iout,'(i5,20f10.5)')
754 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
757 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
758 do l=1,ntermd_2(i,j,k,iblock)
759 write (iout,'(i5,20f10.5)')
760 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
761 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
770 C Read of Side-chain backbone correlation parameters
771 C Modified 11 May 2012 by Adasko
774 read (isccor,*,end=119,err=119) nsccortyp
776 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
778 isccortyp(i)=-isccortyp(-i)
780 iscprol=isccortyp(20)
781 c write (iout,*) 'ntortyp',ntortyp
783 cc maxinter is maximum interaction sites
787 read (isccor,*,end=119,err=119)
788 &nterm_sccor(i,j),nlor_sccor(i,j)
794 nterm_sccor(-i,j)=nterm_sccor(i,j)
795 nterm_sccor(-i,-j)=nterm_sccor(i,j)
796 nterm_sccor(i,-j)=nterm_sccor(i,j)
797 do k=1,nterm_sccor(i,j)
798 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
800 if (j.eq.iscprol) then
801 if (i.eq.isccortyp(10)) then
802 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
803 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
805 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
806 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
807 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
808 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
809 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
810 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
811 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
812 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
815 if (i.eq.isccortyp(10)) then
816 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
817 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
819 if (j.eq.isccortyp(10)) then
820 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
821 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
823 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
824 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
825 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
826 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
827 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
828 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
832 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
833 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
834 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
835 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
838 do k=1,nlor_sccor(i,j)
839 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
840 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
841 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
842 &(1+vlor3sccor(k,i,j)**2)
844 v0sccor(l,i,j)=v0ijsccor
845 v0sccor(l,-i,j)=v0ijsccor1
846 v0sccor(l,i,-j)=v0ijsccor2
847 v0sccor(l,-i,-j)=v0ijsccor3
853 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
854 c write (iout,*) 'ntortyp',ntortyp
856 cc maxinter is maximum interaction sites
860 read (isccor,*,end=119,err=119)
861 & nterm_sccor(i,j),nlor_sccor(i,j)
865 do k=1,nterm_sccor(i,j)
866 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
868 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
871 do k=1,nlor_sccor(i,j)
872 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
873 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
874 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
875 &(1+vlor3sccor(k,i,j)**2)
877 v0sccor(i,j,iblock)=v0ijsccor
885 write (iout,'(/a/)') 'Torsional constants:'
888 write (iout,*) 'ityp',i,' jtyp',j
889 write (iout,*) 'Fourier constants'
890 do k=1,nterm_sccor(i,j)
891 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
893 write (iout,*) 'Lorenz constants'
894 do k=1,nlor_sccor(i,j)
895 write (iout,'(3(1pe15.5))')
896 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
903 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
904 C interaction energy of the Gly, Ala, and Pro prototypes.
908 write (iout,*) "Coefficients of the cumulants"
910 read (ifourier,*) nloctyp
912 read (ifourier,*,end=115,err=115)
913 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
915 write (iout,*) 'Type',i
916 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
958 c Ctilde(1,1,i)=0.0d0
959 c Ctilde(1,2,i)=0.0d0
960 c Ctilde(2,1,i)=0.0d0
961 c Ctilde(2,2,i)=0.0d0
983 c Dtilde(1,1,i)=0.0d0
984 c Dtilde(1,2,i)=0.0d0
985 c Dtilde(2,1,i)=0.0d0
986 c Dtilde(2,2,i)=0.0d0
987 EE(1,1,i)= b(10)+b(11)
988 EE(2,2,i)=-b(10)+b(11)
989 EE(2,1,i)= b(12)-b(13)
990 EE(1,2,i)= b(12)+b(13)
991 EE(1,1,-i)= b(10)+b(11)
992 EE(2,2,-i)=-b(10)+b(11)
993 EE(2,1,-i)=-b(12)+b(13)
994 EE(1,2,-i)=-b(12)-b(13)
1000 c ee(2,1,i)=ee(1,2,i)
1004 write (iout,*) 'Type',i
1006 write(iout,*) B1(1,i),B1(2,i)
1008 write(iout,*) B2(1,i),B2(2,i)
1011 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1015 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1019 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1025 C Read electrostatic-interaction parameters
1029 write (iout,'(/a)') 'Electrostatic interaction constants:'
1030 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1031 & 'IT','JT','APP','BPP','AEL6','AEL3'
1033 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1034 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1035 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1036 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1041 app (i,j)=epp(i,j)*rri*rri
1042 bpp (i,j)=-2.0D0*epp(i,j)*rri
1043 ael6(i,j)=elpp6(i,j)*4.2D0**6
1044 ael3(i,j)=elpp3(i,j)*4.2D0**3
1046 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1047 & ael6(i,j),ael3(i,j)
1052 C Read side-chain interaction parameters.
1054 read (isidep,*,end=117,err=117) ipot,expon
1055 if (ipot.lt.1 .or. ipot.gt.5) then
1056 write (iout,'(2a)') 'Error while reading SC interaction',
1057 & 'potential file - unknown potential type.'
1059 call MPI_Finalize(Ierror)
1065 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1066 & ', exponents are ',expon,2*expon
1067 goto (10,20,30,30,40) ipot
1068 C----------------------- LJ potential ---------------------------------
1069 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1070 & (sigma0(i),i=1,ntyp)
1072 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1073 write (iout,'(a/)') 'The epsilon array:'
1074 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1075 write (iout,'(/a)') 'One-body parameters:'
1076 write (iout,'(a,4x,a)') 'residue','sigma'
1077 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1080 C----------------------- LJK potential --------------------------------
1081 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1082 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1084 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1085 write (iout,'(a/)') 'The epsilon array:'
1086 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1087 write (iout,'(/a)') 'One-body parameters:'
1088 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1089 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1093 C---------------------- GB or BP potential -----------------------------
1095 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1097 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
1098 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
1099 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
1100 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
1101 C For the GB potential convert sigma'**2 into chi'
1104 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1108 write (iout,'(/a/)') 'Parameters of the BP potential:'
1109 write (iout,'(a/)') 'The epsilon array:'
1110 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1111 write (iout,'(/a)') 'One-body parameters:'
1112 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1114 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1115 & chip(i),alp(i),i=1,ntyp)
1118 C--------------------- GBV potential -----------------------------------
1119 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1120 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1121 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1123 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1124 write (iout,'(a/)') 'The epsilon array:'
1125 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1126 write (iout,'(/a)') 'One-body parameters:'
1127 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1128 & 's||/s_|_^2',' chip ',' alph '
1129 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1130 & sigii(i),chip(i),alp(i),i=1,ntyp)
1134 C-----------------------------------------------------------------------
1135 C Calculate the "working" parameters of SC interactions.
1143 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1144 sigma(j,i)=sigma(i,j)
1145 rs0(i,j)=dwa16*sigma(i,j)
1149 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1150 & 'Working parameters of the SC interactions:',
1151 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1156 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1165 sigeps=dsign(1.0D0,epsij)
1167 aa(i,j)=epsij*rrij*rrij
1168 bb(i,j)=-sigeps*epsij*rrij
1172 sigt1sq=sigma0(i)**2
1173 sigt2sq=sigma0(j)**2
1176 ratsig1=sigt2sq/sigt1sq
1177 ratsig2=1.0D0/ratsig1
1178 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1179 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1180 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1184 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1185 sigmaii(i,j)=rsum_max
1186 sigmaii(j,i)=rsum_max
1188 c sigmaii(i,j)=r0(i,j)
1189 c sigmaii(j,i)=r0(i,j)
1191 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1192 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1193 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1194 augm(i,j)=epsij*r_augm**(2*expon)
1195 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1202 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1203 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1204 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1210 C Define the SC-p interaction constants (hard-coded; old style)
1213 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1215 c aad(i,1)=0.3D0*4.0D0**12
1216 C Following line for constants currently implemented
1217 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1218 aad(i,1)=1.5D0*4.0D0**12
1219 c aad(i,1)=0.17D0*5.6D0**12
1221 C "Soft" SC-p repulsion
1223 C Following line for constants currently implemented
1224 c aad(i,1)=0.3D0*4.0D0**6
1225 C "Hard" SC-p repulsion
1226 bad(i,1)=3.0D0*4.0D0**6
1227 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1236 C 8/9/01 Read the SC-p interaction constants from file
1239 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1242 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1243 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1244 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1245 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1249 write (iout,*) "Parameters of SC-p interactions:"
1251 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1252 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1258 C Define the constants of the disulfide bridge
1262 c Old arbitrary potential - commented out.
1267 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1268 c energy surface of diethyl disulfide.
1269 c A. Liwo and U. Kozlowska, 11/24/03
1285 C if(me.eq.king) then
1286 C write (iout,'(/a)') "Disulfide bridge parameters:"
1287 C write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1288 C write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1289 C write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1290 C write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1294 111 write (iout,*) "Error reading bending energy parameters."
1296 112 write (iout,*) "Error reading rotamer energy parameters."
1298 113 write (iout,*) "Error reading torsional energy parameters."
1300 114 write (iout,*) "Error reading double torsional energy parameters."
1303 & "Error reading cumulant (multibody energy) parameters."
1305 116 write (iout,*) "Error reading electrostatic energy parameters."
1307 117 write (iout,*) "Error reading side chain interaction parameters."
1309 118 write (iout,*) "Error reading SCp interaction parameters."
1311 119 write (iout,*) "Error reading SCCOR parameters"
1314 call MPI_Finalize(Ierror)
1321 subroutine getenv_loc(var, val)
1322 character(*) var, val
1325 character(2000) line
1328 open (196,file='env',status='old',readonly,shared)
1330 c write(*,*)'looking for ',var
1331 10 read(196,*,err=11,end=11)line
1332 iread=index(line,var)
1333 c write(*,*)iread,' ',var,' ',line
1334 if (iread.eq.0) go to 10
1335 c write(*,*)'---> ',line
1341 iread=iread+ilen(var)+1
1342 read (line(iread:),*,err=12,end=12) val
1343 c write(*,*)'OK: ',var,' = ',val
1349 #elif (defined CRAY)
1350 integer lennam,lenval,ierror
1352 c getenv using a POSIX call, useful on the T3D
1353 c Sept 1996, comment out error check on advice of H. Pritchard
1356 if(lennam.le.0) stop '--error calling getenv--'
1357 call pxfgetenv(var,lennam,val,lenval,ierror)
1358 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1360 call getenv(var,val)