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 write (iout,*) "TUTUTU"
916 read (ifourier,*,end=115,err=115) (bnew1(ii,1,i),ii=1,3)
917 read (ifourier,*,end=115,err=115) (bnew2(ii,1,i),ii=1,3)
918 read (ifourier,*,end=115,err=115) (bnew1(ii,2,i),ii=1,1)
919 read (ifourier,*,end=115,err=115) (bnew2(ii,2,i),ii=1,1)
920 read (ifourier,*,end=115,err=115) (eenew(ii,i),ii=1,1)
924 write (iout,*) 'Type',i
925 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
933 c B1tilde(1,i) = b(3)
934 c B1tilde(2,i) =-b(5)
935 c B1tilde(1,-i) =-b(3)
936 c B1tilde(2,-i) =b(5)
967 c Ctilde(1,1,i)=0.0d0
968 c Ctilde(1,2,i)=0.0d0
969 c Ctilde(2,1,i)=0.0d0
970 c Ctilde(2,2,i)=0.0d0
992 c Dtilde(1,1,i)=0.0d0
993 c Dtilde(1,2,i)=0.0d0
994 c Dtilde(2,1,i)=0.0d0
995 c Dtilde(2,2,i)=0.0d0
996 EEold(1,1,i)= b(10)+b(11)
997 EEold(2,2,i)=-b(10)+b(11)
998 EEold(2,1,i)= b(12)-b(13)
999 EEold(1,2,i)= b(12)+b(13)
1000 EEold(1,1,-i)= b(10)+b(11)
1001 EEold(2,2,-i)=-b(10)+b(11)
1002 EEold(2,1,-i)=-b(12)+b(13)
1003 EEold(1,2,-i)=-b(12)-b(13)
1004 write(iout,*) "TU DOCHODZE"
1009 c ee(2,1,i)=ee(1,2,i)
1014 write (iout,*) 'Type',i
1016 write(iout,*) B1(1,i),B1(2,i)
1018 write(iout,*) B2(1,i),B2(2,i)
1021 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1025 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1029 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1036 C Read electrostatic-interaction parameters
1040 write (iout,'(/a)') 'Electrostatic interaction constants:'
1041 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1042 & 'IT','JT','APP','BPP','AEL6','AEL3'
1044 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1045 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1046 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1047 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1052 app (i,j)=epp(i,j)*rri*rri
1053 bpp (i,j)=-2.0D0*epp(i,j)*rri
1054 ael6(i,j)=elpp6(i,j)*4.2D0**6
1055 ael3(i,j)=elpp3(i,j)*4.2D0**3
1057 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1058 & ael6(i,j),ael3(i,j)
1063 C Read side-chain interaction parameters.
1065 read (isidep,*,end=117,err=117) ipot,expon
1066 if (ipot.lt.1 .or. ipot.gt.5) then
1067 write (iout,'(2a)') 'Error while reading SC interaction',
1068 & 'potential file - unknown potential type.'
1070 call MPI_Finalize(Ierror)
1076 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1077 & ', exponents are ',expon,2*expon
1078 goto (10,20,30,30,40) ipot
1079 C----------------------- LJ potential ---------------------------------
1080 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1081 & (sigma0(i),i=1,ntyp)
1083 write (iout,'(/a/)') 'Parameters of the LJ 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,a)') 'residue','sigma'
1088 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1091 C----------------------- LJK potential --------------------------------
1092 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1093 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1095 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1096 write (iout,'(a/)') 'The epsilon array:'
1097 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1098 write (iout,'(/a)') 'One-body parameters:'
1099 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1100 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1104 C---------------------- GB or BP potential -----------------------------
1106 read (isidep,*,end=116,err=116)(eps(i,j),j=i,ntyp)
1108 read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
1109 read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
1110 read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
1111 read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
1112 C For the GB potential convert sigma'**2 into chi'
1115 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1119 write (iout,'(/a/)') 'Parameters of the BP potential:'
1120 write (iout,'(a/)') 'The epsilon array:'
1121 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1122 write (iout,'(/a)') 'One-body parameters:'
1123 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1125 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1126 & chip(i),alp(i),i=1,ntyp)
1129 C--------------------- GBV potential -----------------------------------
1130 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1131 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1132 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1134 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1135 write (iout,'(a/)') 'The epsilon array:'
1136 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1137 write (iout,'(/a)') 'One-body parameters:'
1138 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1139 & 's||/s_|_^2',' chip ',' alph '
1140 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1141 & sigii(i),chip(i),alp(i),i=1,ntyp)
1145 C-----------------------------------------------------------------------
1146 C Calculate the "working" parameters of SC interactions.
1154 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1155 sigma(j,i)=sigma(i,j)
1156 rs0(i,j)=dwa16*sigma(i,j)
1160 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1161 & 'Working parameters of the SC interactions:',
1162 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1167 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1176 sigeps=dsign(1.0D0,epsij)
1178 aa(i,j)=epsij*rrij*rrij
1179 bb(i,j)=-sigeps*epsij*rrij
1183 sigt1sq=sigma0(i)**2
1184 sigt2sq=sigma0(j)**2
1187 ratsig1=sigt2sq/sigt1sq
1188 ratsig2=1.0D0/ratsig1
1189 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1190 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1191 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1195 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1196 sigmaii(i,j)=rsum_max
1197 sigmaii(j,i)=rsum_max
1199 c sigmaii(i,j)=r0(i,j)
1200 c sigmaii(j,i)=r0(i,j)
1202 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1203 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1204 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1205 augm(i,j)=epsij*r_augm**(2*expon)
1206 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1213 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1214 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1215 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1221 C Define the SC-p interaction constants (hard-coded; old style)
1224 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1226 c aad(i,1)=0.3D0*4.0D0**12
1227 C Following line for constants currently implemented
1228 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1229 aad(i,1)=1.5D0*4.0D0**12
1230 c aad(i,1)=0.17D0*5.6D0**12
1232 C "Soft" SC-p repulsion
1234 C Following line for constants currently implemented
1235 c aad(i,1)=0.3D0*4.0D0**6
1236 C "Hard" SC-p repulsion
1237 bad(i,1)=3.0D0*4.0D0**6
1238 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1247 C 8/9/01 Read the SC-p interaction constants from file
1250 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1253 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1254 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1255 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1256 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1260 write (iout,*) "Parameters of SC-p interactions:"
1262 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1263 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1269 C Define the constants of the disulfide bridge
1273 c Old arbitrary potential - commented out.
1278 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1279 c energy surface of diethyl disulfide.
1280 c A. Liwo and U. Kozlowska, 11/24/03
1297 write (iout,'(/a)') "Disulfide bridge parameters:"
1298 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1299 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1300 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1301 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1305 111 write (iout,*) "Error reading bending energy parameters."
1307 112 write (iout,*) "Error reading rotamer energy parameters."
1309 113 write (iout,*) "Error reading torsional energy parameters."
1311 114 write (iout,*) "Error reading double torsional energy parameters."
1314 & "Error reading cumulant (multibody energy) parameters."
1316 116 write (iout,*) "Error reading electrostatic energy parameters."
1318 117 write (iout,*) "Error reading side chain interaction parameters."
1320 118 write (iout,*) "Error reading SCp interaction parameters."
1322 119 write (iout,*) "Error reading SCCOR parameters"
1325 call MPI_Finalize(Ierror)
1332 subroutine getenv_loc(var, val)
1333 character(*) var, val
1336 character(2000) line
1339 open (196,file='env',status='old',readonly,shared)
1341 c write(*,*)'looking for ',var
1342 10 read(196,*,err=11,end=11)line
1343 iread=index(line,var)
1344 c write(*,*)iread,' ',var,' ',line
1345 if (iread.eq.0) go to 10
1346 c write(*,*)'---> ',line
1352 iread=iread+ilen(var)+1
1353 read (line(iread:),*,err=12,end=12) val
1354 c write(*,*)'OK: ',var,' = ',val
1360 #elif (defined CRAY)
1361 integer lennam,lenval,ierror
1363 c getenv using a POSIX call, useful on the T3D
1364 c Sept 1996, comment out error check on advice of H. Pritchard
1367 if(lennam.le.0) stop '--error calling getenv--'
1368 call pxfgetenv(var,lennam,val,lenval,ierror)
1369 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1371 call getenv(var,val)