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 EEold(1,1,i)= b(10)+b(11)
995 EEold(2,2,i)=-b(10)+b(11)
996 EEold(2,1,i)= b(12)-b(13)
997 EEold(1,2,i)= b(12)+b(13)
998 EEold(1,1,-i)= b(10)+b(11)
999 EEold(2,2,-i)=-b(10)+b(11)
1000 EEold(2,1,-i)=-b(12)+b(13)
1001 EEold(1,2,-i)=-b(12)-b(13)
1002 write(iout,*) "TU DOCHODZE"
1007 c ee(2,1,i)=ee(1,2,i)
1012 write (iout,*) 'Type',i
1014 write(iout,*) B1(1,i),B1(2,i)
1016 write(iout,*) B2(1,i),B2(2,i)
1019 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1023 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1027 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1034 C Read electrostatic-interaction parameters
1038 write (iout,'(/a)') 'Electrostatic interaction constants:'
1039 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1040 & 'IT','JT','APP','BPP','AEL6','AEL3'
1042 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1043 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1044 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1045 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1050 app (i,j)=epp(i,j)*rri*rri
1051 bpp (i,j)=-2.0D0*epp(i,j)*rri
1052 ael6(i,j)=elpp6(i,j)*4.2D0**6
1053 ael3(i,j)=elpp3(i,j)*4.2D0**3
1055 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1056 & ael6(i,j),ael3(i,j)
1061 C Read side-chain interaction parameters.
1063 read (isidep,*,end=117,err=117) ipot,expon
1064 if (ipot.lt.1 .or. ipot.gt.5) then
1065 write (iout,'(2a)') 'Error while reading SC interaction',
1066 & 'potential file - unknown potential type.'
1068 call MPI_Finalize(Ierror)
1074 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1075 & ', exponents are ',expon,2*expon
1076 goto (10,20,30,30,40) ipot
1077 C----------------------- LJ potential ---------------------------------
1078 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1079 & (sigma0(i),i=1,ntyp)
1081 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1082 write (iout,'(a/)') 'The epsilon array:'
1083 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1084 write (iout,'(/a)') 'One-body parameters:'
1085 write (iout,'(a,4x,a)') 'residue','sigma'
1086 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1089 C----------------------- LJK potential --------------------------------
1090 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1091 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1093 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1094 write (iout,'(a/)') 'The epsilon array:'
1095 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1096 write (iout,'(/a)') 'One-body parameters:'
1097 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1098 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1102 C---------------------- GB or BP potential -----------------------------
1104 read (isidep,*,end=116,err=116)(eps(i,j),j=i,ntyp)
1106 read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
1107 read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
1108 read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
1109 read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
1110 C For the GB potential convert sigma'**2 into chi'
1113 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1117 write (iout,'(/a/)') 'Parameters of the BP potential:'
1118 write (iout,'(a/)') 'The epsilon array:'
1119 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1120 write (iout,'(/a)') 'One-body parameters:'
1121 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1123 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1124 & chip(i),alp(i),i=1,ntyp)
1127 C--------------------- GBV potential -----------------------------------
1128 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1129 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1130 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1132 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1133 write (iout,'(a/)') 'The epsilon array:'
1134 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1135 write (iout,'(/a)') 'One-body parameters:'
1136 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1137 & 's||/s_|_^2',' chip ',' alph '
1138 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1139 & sigii(i),chip(i),alp(i),i=1,ntyp)
1143 C-----------------------------------------------------------------------
1144 C Calculate the "working" parameters of SC interactions.
1152 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1153 sigma(j,i)=sigma(i,j)
1154 rs0(i,j)=dwa16*sigma(i,j)
1158 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1159 & 'Working parameters of the SC interactions:',
1160 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1165 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1174 sigeps=dsign(1.0D0,epsij)
1176 aa(i,j)=epsij*rrij*rrij
1177 bb(i,j)=-sigeps*epsij*rrij
1181 sigt1sq=sigma0(i)**2
1182 sigt2sq=sigma0(j)**2
1185 ratsig1=sigt2sq/sigt1sq
1186 ratsig2=1.0D0/ratsig1
1187 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1188 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1189 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1193 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1194 sigmaii(i,j)=rsum_max
1195 sigmaii(j,i)=rsum_max
1197 c sigmaii(i,j)=r0(i,j)
1198 c sigmaii(j,i)=r0(i,j)
1200 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1201 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1202 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1203 augm(i,j)=epsij*r_augm**(2*expon)
1204 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1211 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1212 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1213 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1219 C Define the SC-p interaction constants (hard-coded; old style)
1222 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1224 c aad(i,1)=0.3D0*4.0D0**12
1225 C Following line for constants currently implemented
1226 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1227 aad(i,1)=1.5D0*4.0D0**12
1228 c aad(i,1)=0.17D0*5.6D0**12
1230 C "Soft" SC-p repulsion
1232 C Following line for constants currently implemented
1233 c aad(i,1)=0.3D0*4.0D0**6
1234 C "Hard" SC-p repulsion
1235 bad(i,1)=3.0D0*4.0D0**6
1236 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1245 C 8/9/01 Read the SC-p interaction constants from file
1248 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1251 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1252 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1253 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1254 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1258 write (iout,*) "Parameters of SC-p interactions:"
1260 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1261 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1267 C Define the constants of the disulfide bridge
1271 c Old arbitrary potential - commented out.
1276 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1277 c energy surface of diethyl disulfide.
1278 c A. Liwo and U. Kozlowska, 11/24/03
1295 write (iout,'(/a)') "Disulfide bridge parameters:"
1296 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1297 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1298 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1299 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1303 111 write (iout,*) "Error reading bending energy parameters."
1305 112 write (iout,*) "Error reading rotamer energy parameters."
1307 113 write (iout,*) "Error reading torsional energy parameters."
1309 114 write (iout,*) "Error reading double torsional energy parameters."
1312 & "Error reading cumulant (multibody energy) parameters."
1314 116 write (iout,*) "Error reading electrostatic energy parameters."
1316 117 write (iout,*) "Error reading side chain interaction parameters."
1318 118 write (iout,*) "Error reading SCp interaction parameters."
1320 119 write (iout,*) "Error reading SCCOR parameters"
1323 call MPI_Finalize(Ierror)
1330 subroutine getenv_loc(var, val)
1331 character(*) var, val
1334 character(2000) line
1337 open (196,file='env',status='old',readonly,shared)
1339 c write(*,*)'looking for ',var
1340 10 read(196,*,err=11,end=11)line
1341 iread=index(line,var)
1342 c write(*,*)iread,' ',var,' ',line
1343 if (iread.eq.0) go to 10
1344 c write(*,*)'---> ',line
1350 iread=iread+ilen(var)+1
1351 read (line(iread:),*,err=12,end=12) val
1352 c write(*,*)'OK: ',var,' = ',val
1358 #elif (defined CRAY)
1359 integer lennam,lenval,ierror
1361 c getenv using a POSIX call, useful on the T3D
1362 c Sept 1996, comment out error check on advice of H. Pritchard
1365 if(lennam.le.0) stop '--error calling getenv--'
1366 call pxfgetenv(var,lennam,val,lenval,ierror)
1367 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1369 call getenv(var,val)