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 C write(iout,*) "I am here",ntyp1
215 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
216 C write(iout,*) "I am herew"
218 ithetyp(i)=-ithetyp(-i)
221 do i=-maxthetyp,maxthetyp
222 do j=-maxthetyp,maxthetyp
223 do k=-maxthetyp,maxthetyp
224 aa0thet(i,j,k,iblock)=0.0d0
226 aathet(l,i,j,k,iblock)=0.0d0
230 bbthet(m,l,i,j,k,iblock)=0.0d0
231 ccthet(m,l,i,j,k,iblock)=0.0d0
232 ddthet(m,l,i,j,k,iblock)=0.0d0
233 eethet(m,l,i,j,k,iblock)=0.0d0
239 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
240 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
248 c VAR:iblock means terminally blocking group 1=non-proline 2=proline
250 c VAR:ntethtyp is type of theta potentials type currently 0=glycine
251 c VAR:1=non-glicyne non-proline 2=proline
252 c VAR:negative values for D-aminoacid
254 do j=-nthetyp,nthetyp
255 do k=-nthetyp,nthetyp
256 read (ithep,'(6a)',end=111,err=111) res1
257 c VAR: aa0thet is variable describing the average value of Foureir
258 c VAR: expansion series
259 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
260 c VAR: aathet is foureir expansion in theta/2 angle for full formula
261 c VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
262 Condens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
263 read (ithep,*,end=111,err=111)
264 &(aathet(l,i,j,k,iblock),l=1,ntheterm)
265 read (ithep,*,end=111,err=111)
266 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
267 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
268 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
269 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),
271 read (ithep,*,end=111,err=111)
272 & (((ffthet(llll,lll,ll,i,j,k,iblock),
273 & ffthet(lll,llll,ll,i,j,k,iblock),
274 & ggthet(llll,lll,ll,i,j,k,iblock),
275 & ggthet(lll,llll,ll,i,j,k,iblock),
276 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
283 C For dummy ends assign glycine-type coefficients of theta-only terms; the
284 C coefficients of theta-and-gamma-dependent terms are zero.
285 C IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
286 C RECOMENTDED AFTER VERSION 3.3)
290 c aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
291 c aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
293 c aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
294 c aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
297 c aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
299 c aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
302 C AND COMMENT THE LOOPS BELOW
306 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
307 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
309 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
310 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
313 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
315 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
319 C Substitution for D aminoacids from symmetry.
322 do j=-nthetyp,nthetyp
323 do k=-nthetyp,nthetyp
324 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
326 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
330 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
331 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
332 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
333 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
339 ffthet(llll,lll,ll,i,j,k,iblock)=
340 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
341 ffthet(lll,llll,ll,i,j,k,iblock)=
342 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
343 ggthet(llll,lll,ll,i,j,k,iblock)=
344 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
345 ggthet(lll,llll,ll,i,j,k,iblock)=
346 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
355 C Control printout of the coefficients of virtual-bond-angle potentials
358 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
362 write (iout,'(//4a)')
363 & 'Type ',onelett(i),onelett(j),onelett(k)
364 write (iout,'(//a,10x,a)') " l","a[l]"
365 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
366 write (iout,'(i2,1pe15.5)')
367 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
369 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))')
370 & "b",l,"c",l,"d",l,"e",l
372 write (iout,'(i2,4(1pe15.5))') m,
373 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
374 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
378 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
379 & "f+",l,"f-",l,"g+",l,"g-",l
382 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
383 & ffthet(n,m,l,i,j,k,iblock),
384 & ffthet(m,n,l,i,j,k,iblock),
385 & ggthet(n,m,l,i,j,k,iblock),
386 & ggthet(m,n,l,i,j,k,iblock)
395 write (2,*) "Start reading THETA_PDB"
397 read (ithep_pdb,*,err=111,end=111)
398 & a0thet(i),(athet(j,i,1,1),j=1,2),
399 & (bthet(j,i,1,1),j=1,2)
400 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
401 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
402 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
406 athet(1,i,1,-1)=athet(1,i,1,1)
407 athet(2,i,1,-1)=athet(2,i,1,1)
408 bthet(1,i,1,-1)=-bthet(1,i,1,1)
409 bthet(2,i,1,-1)=-bthet(2,i,1,1)
410 athet(1,i,-1,1)=-athet(1,i,1,1)
411 athet(2,i,-1,1)=-athet(2,i,1,1)
412 bthet(1,i,-1,1)=bthet(1,i,1,1)
413 bthet(2,i,-1,1)=bthet(2,i,1,1)
417 athet(1,i,-1,-1)=athet(1,-i,1,1)
418 athet(2,i,-1,-1)=-athet(2,-i,1,1)
419 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
420 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
421 athet(1,i,-1,1)=athet(1,-i,1,1)
422 athet(2,i,-1,1)=-athet(2,-i,1,1)
423 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
424 bthet(2,i,-1,1)=bthet(2,-i,1,1)
425 athet(1,i,1,-1)=-athet(1,-i,1,1)
426 athet(2,i,1,-1)=athet(2,-i,1,1)
427 bthet(1,i,1,-1)=bthet(1,-i,1,1)
428 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
433 polthet(j,i)=polthet(j,-i)
436 gthet(j,i)=gthet(j,-i)
439 write (2,*) "End reading THETA_PDB"
445 C Read the parameters of the probability distribution/energy expression
446 C of the side chains.
449 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
453 dsc_inv(i)=1.0D0/dsc(i)
464 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
465 & ((blower(k,l,1),l=1,k),k=1,3)
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)
476 akl=akl+blower(k,m,j)*blower(l,m,j)
487 write (iout,'(/a)') 'Parameters of side-chain local geometry'
492 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
493 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
494 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
495 & 'log h',(bsc(j,i),j=1,nlobi)
496 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
497 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
499 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
500 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
503 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
504 write (iout,'(a,f10.4,4(16x,f10.4))')
505 & 'Center ',(bsc(j,i),j=1,nlobi)
506 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
515 C Read scrot parameters for potentials determined from all-atom AM1 calculations
516 C added by Urszula Kozlowska 07/11/2007
519 read (irotam,*,end=112,err=112)
521 read (irotam,*,end=112,err=112)
524 read(irotam,*,end=112,err=112) sc_parmin(j,i)
529 C Read the parameters of the probability distribution/energy expression
530 C of the side chains.
532 write (2,*) "Start reading ROTAM_PDB"
535 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
539 dsc_inv(i)=1.0D0/dsc(i)
550 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
551 & ((blower(k,l,1),l=1,k),k=1,3)
553 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
554 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
555 & ((blower(k,l,j),l=1,k),k=1,3)
562 akl=akl+blower(k,m,j)*blower(l,m,j)
572 write (2,*) "Ending reading ROTAM_PDB"
578 C Read torsional parameters in old format
580 read (itorp,*,end=113,err=113) ntortyp,nterm_old
581 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
582 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
587 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
593 write (iout,'(/a/)') 'Torsional constants:'
596 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
597 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
603 C Read torsional parameters
605 read (itorp,*,end=113,err=113) ntortyp
606 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
609 itortyp(i)=-itortyp(-i)
611 write (iout,*) 'ntortyp',ntortyp
613 do j=-ntortyp+1,ntortyp-1
614 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
616 nterm(-i,-j,iblock)=nterm(i,j,iblock)
617 nlor(-i,-j,iblock)=nlor(i,j,iblock)
620 do k=1,nterm(i,j,iblock)
621 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
623 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
624 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
625 v0ij=v0ij+si*v1(k,i,j,iblock)
627 c write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
628 c write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
629 c &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
631 do k=1,nlor(i,j,iblock)
632 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
633 & vlor2(k,i,j),vlor3(k,i,j)
634 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
637 v0(-i,-j,iblock)=v0ij
643 write (iout,'(/a/)') 'Torsional constants:'
646 write (iout,*) 'ityp',i,' jtyp',j
647 write (iout,*) 'Fourier constants'
648 do k=1,nterm(i,j,iblock)
649 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
652 write (iout,*) 'Lorenz constants'
653 do k=1,nlor(i,j,iblock)
654 write (iout,'(3(1pe15.5))')
655 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
661 C 6/23/01 Read parameters for double torsionals
668 do j=-ntortyp+1,ntortyp-1
669 do k=-ntortyp+1,ntortyp-1
670 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
671 c write (iout,*) "OK onelett",
674 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
675 & .or. t3.ne.toronelet(k)) then
676 write (iout,*) "Error in double torsional parameter file",
679 call MPI_Finalize(Ierror)
681 stop "Error in double torsional parameter file"
683 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
684 & ntermd_2(i,j,k,iblock)
685 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
686 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
687 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
688 & ntermd_1(i,j,k,iblock))
689 read (itordp,*,end=114,err=114)(v1s(1,l,i,j,k,iblock),l=1,
690 & ntermd_1(i,j,k,iblock))
691 read (itordp,*,end=114,err=114)(v1c(2,l,i,j,k,iblock),l=1,
692 & ntermd_1(i,j,k,iblock))
693 read (itordp,*,end=114,err=114)(v1s(2,l,i,j,k,iblock),l=1,
694 & ntermd_1(i,j,k,iblock))
695 C Martix of D parameters for one dimesional foureir series
696 do l=1,ntermd_1(i,j,k,iblock)
697 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
698 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
699 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
700 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
701 c write(iout,*) "whcodze" ,
702 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
704 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
705 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
706 & v2s(m,l,i,j,k,iblock),
707 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
708 C Martix of D parameters for two dimesional fourier series
709 do l=1,ntermd_2(i,j,k,iblock)
711 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
712 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
713 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
714 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
723 write (iout,*) 'Constants for double torsionals'
726 do j=-ntortyp+1,ntortyp-1
727 do k=-ntortyp+1,ntortyp-1
728 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
729 & ' nsingle',ntermd_1(i,j,k,iblock),
730 & ' ndouble',ntermd_2(i,j,k,iblock)
732 write (iout,*) 'Single angles:'
733 do l=1,ntermd_1(i,j,k,iblock)
734 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
735 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
736 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
737 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
740 write (iout,*) 'Pairs of angles:'
741 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
742 do l=1,ntermd_2(i,j,k,iblock)
743 write (iout,'(i5,20f10.5)')
744 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
747 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
748 do l=1,ntermd_2(i,j,k,iblock)
749 write (iout,'(i5,20f10.5)')
750 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
751 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
760 C Read of Side-chain backbone correlation parameters
761 C Modified 11 May 2012 by Adasko
764 read (isccor,*,end=1113,err=1113) nsccortyp
766 write (iout,*) "Tu wchodze"
767 read (isccor,*,end=1113,err=1113) (isccortyp(i),i=1,ntyp)
769 isccortyp(i)=-isccortyp(-i)
771 iscprol=isccortyp(20)
772 C write (iout,*) 'ntortyp',ntortyp
774 cc maxinter is maximum interaction sites
778 read (isccor,*,end=1113,err=1113) nterm_sccor(i,j),
785 nterm_sccor(-i,j)=nterm_sccor(i,j)
786 nterm_sccor(-i,-j)=nterm_sccor(i,j)
787 nterm_sccor(i,-j)=nterm_sccor(i,j)
788 do k=1,nterm_sccor(i,j)
789 read (isccor,*,end=1113,err=1113) kk,v1sccor(k,l,i,j)
791 c write(iout,*) "k=",kk
792 if (j.eq.iscprol) then
793 if (i.eq.isccortyp(10)) then
794 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
795 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
797 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
798 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
799 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
800 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
801 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
802 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
803 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
804 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
807 if (i.eq.isccortyp(10)) then
808 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
809 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
811 if (j.eq.isccortyp(10)) then
812 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
813 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
815 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
816 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
817 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
818 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
819 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
820 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
824 C read (isccor,*,end=1113,err=1113) kk,v1sccor(k,l,i,j)
825 C & ,v2sccor(k,l,i,j)
826 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
827 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
828 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
829 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
832 do k=1,nlor_sccor(i,j)
833 read (isccor,*,end=1113,err=1113) kk,vlor1sccor(k,i,j),
834 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
835 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
836 &(1+vlor3sccor(k,i,j)**2)
838 v0sccor(l,i,j)=v0ijsccor
839 v0sccor(l,-i,j)=v0ijsccor1
840 v0sccor(l,i,-j)=v0ijsccor2
841 v0sccor(l,-i,-j)=v0ijsccor3
847 write(iout,*) "a tu nie wchodze"
848 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
849 c write (iout,*) 'ntortyp',ntortyp
851 cc maxinter is maximum interaction sites
855 read (isccor,*,end=119,err=119)
856 & nterm_sccor(i,j),nlor_sccor(i,j)
860 do k=1,nterm_sccor(i,j)
861 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
863 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
866 do k=1,nlor_sccor(i,j)
867 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
868 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
869 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
870 &(1+vlor3sccor(k,i,j)**2)
872 v0sccor(i,j,iblock)=v0ijsccor
880 write (iout,'(/a/)') 'Torsional constants:'
883 write (iout,*) 'ityp',i,' jtyp',j
884 write (iout,*) 'Fourier constants'
885 do k=1,nterm_sccor(i,j)
886 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
888 write (iout,*) 'Lorenz constants'
889 do k=1,nlor_sccor(i,j)
890 write (iout,'(3(1pe15.5))')
891 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
898 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
899 C interaction energy of the Gly, Ala, and Pro prototypes.
903 write (iout,*) "Coefficients of the cumulants"
905 read (ifourier,*) nloctyp
907 read (ifourier,*,end=115,err=115)
908 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
910 write (iout,*) 'Type',i
911 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
939 c Ctilde(1,1,i)=0.0d0
940 c Ctilde(1,2,i)=0.0d0
941 c Ctilde(2,1,i)=0.0d0
942 c Ctilde(2,2,i)=0.0d0
955 c Dtilde(1,1,i)=0.0d0
956 c Dtilde(1,2,i)=0.0d0
957 c Dtilde(2,1,i)=0.0d0
958 c Dtilde(2,2,i)=0.0d0
959 EE(1,1,i)= b(10)+b(11)
960 EE(2,2,i)=-b(10)+b(11)
961 EE(2,1,i)= b(12)-b(13)
962 EE(1,2,i)= b(12)+b(13)
967 c ee(2,1,i)=ee(1,2,i)
971 write (iout,*) 'Type',i
973 write(iout,*) B1(1,i),B1(2,i)
975 write(iout,*) B2(1,i),B2(2,i)
978 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
982 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
986 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
991 C Read electrostatic-interaction parameters
995 write (iout,'(/a)') 'Electrostatic interaction constants:'
996 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
997 & 'IT','JT','APP','BPP','AEL6','AEL3'
999 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1000 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1001 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1002 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1007 app (i,j)=epp(i,j)*rri*rri
1008 bpp (i,j)=-2.0D0*epp(i,j)*rri
1009 ael6(i,j)=elpp6(i,j)*4.2D0**6
1010 ael3(i,j)=elpp3(i,j)*4.2D0**3
1011 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1012 & ael6(i,j),ael3(i,j)
1016 C Read side-chain interaction parameters.
1018 read (isidep,*,end=117,err=117) ipot,expon
1019 if (ipot.lt.1 .or. ipot.gt.5) then
1020 write (iout,'(2a)') 'Error while reading SC interaction',
1021 & 'potential file - unknown potential type.'
1023 call MPI_Finalize(Ierror)
1029 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1030 & ', exponents are ',expon,2*expon
1031 goto (10,20,30,30,40) ipot
1032 C----------------------- LJ potential ---------------------------------
1033 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1034 & (sigma0(i),i=1,ntyp)
1036 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1037 write (iout,'(a/)') 'The epsilon array:'
1038 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1039 write (iout,'(/a)') 'One-body parameters:'
1040 write (iout,'(a,4x,a)') 'residue','sigma'
1041 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1044 C----------------------- LJK potential --------------------------------
1045 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1046 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1048 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1049 write (iout,'(a/)') 'The epsilon array:'
1050 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1051 write (iout,'(/a)') 'One-body parameters:'
1052 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1053 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1057 C---------------------- GB or BP potential -----------------------------
1059 read (isidep,*,end=116,err=116)(eps(i,j),j=i,ntyp)
1061 read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
1062 read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
1063 read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
1064 read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
1066 c 30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1067 c & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
1068 c & (alp(i),i=1,ntyp)
1069 C For the GB potential convert sigma'**2 into chi'
1072 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1076 write (iout,'(/a/)') 'Parameters of the BP potential:'
1077 write (iout,'(a/)') 'The epsilon array:'
1078 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1079 write (iout,'(/a)') 'One-body parameters:'
1080 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1082 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1083 & chip(i),alp(i),i=1,ntyp)
1086 C--------------------- GBV potential -----------------------------------
1087 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1088 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1089 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1091 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1092 write (iout,'(a/)') 'The epsilon array:'
1093 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1094 write (iout,'(/a)') 'One-body parameters:'
1095 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1096 & 's||/s_|_^2',' chip ',' alph '
1097 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1098 & sigii(i),chip(i),alp(i),i=1,ntyp)
1102 C-----------------------------------------------------------------------
1103 C Calculate the "working" parameters of SC interactions.
1111 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1112 sigma(j,i)=sigma(i,j)
1113 rs0(i,j)=dwa16*sigma(i,j)
1117 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1118 & 'Working parameters of the SC interactions:',
1119 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1124 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1133 sigeps=dsign(1.0D0,epsij)
1135 aa(i,j)=epsij*rrij*rrij
1136 bb(i,j)=-sigeps*epsij*rrij
1140 sigt1sq=sigma0(i)**2
1141 sigt2sq=sigma0(j)**2
1144 ratsig1=sigt2sq/sigt1sq
1145 ratsig2=1.0D0/ratsig1
1146 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1147 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1148 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1152 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1153 sigmaii(i,j)=rsum_max
1154 sigmaii(j,i)=rsum_max
1156 c sigmaii(i,j)=r0(i,j)
1157 c sigmaii(j,i)=r0(i,j)
1159 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1160 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1161 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1162 augm(i,j)=epsij*r_augm**(2*expon)
1163 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1170 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1171 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1172 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1178 C Define the SC-p interaction constants (hard-coded; old style)
1181 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1183 c aad(i,1)=0.3D0*4.0D0**12
1184 C Following line for constants currently implemented
1185 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1186 aad(i,1)=1.5D0*4.0D0**12
1187 c aad(i,1)=0.17D0*5.6D0**12
1189 C "Soft" SC-p repulsion
1191 C Following line for constants currently implemented
1192 c aad(i,1)=0.3D0*4.0D0**6
1193 C "Hard" SC-p repulsion
1194 bad(i,1)=3.0D0*4.0D0**6
1195 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1204 C 8/9/01 Read the SC-p interaction constants from file
1207 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1210 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1211 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1212 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1213 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1217 write (iout,*) "Parameters of SC-p interactions:"
1219 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1220 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1225 C Define the constants of the disulfide bridge
1229 c Old arbitrary potential - commented out.
1234 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1235 c energy surface of diethyl disulfide.
1236 c A. Liwo and U. Kozlowska, 11/24/03
1253 write (iout,'(/a)') "Disulfide bridge parameters:"
1254 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1255 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1256 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1257 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1261 111 write (iout,*) "Error reading bending energy parameters."
1263 112 write (iout,*) "Error reading rotamer energy parameters."
1265 113 write (iout,*) "Error reading torsional energy parameters."
1268 & "Error reading side-chain torsional energy parameters."
1270 114 write (iout,*) "Error reading double torsional energy parameters."
1273 & "Error reading cumulant (multibody energy) parameters."
1275 116 write (iout,*) "Error reading electrostatic energy parameters."
1277 117 write (iout,*) "Error reading side chain interaction parameters."
1279 118 write (iout,*) "Error reading SCp interaction parameters."
1281 119 write (iout,*) "Error reading SCCOR parameters"
1284 call MPI_Finalize(Ierror)
1291 subroutine getenv_loc(var, val)
1292 character(*) var, val
1295 character(2000) line
1298 open (196,file='env',status='old',readonly,shared)
1300 c write(*,*)'looking for ',var
1301 10 read(196,*,err=11,end=11)line
1302 iread=index(line,var)
1303 c write(*,*)iread,' ',var,' ',line
1304 if (iread.eq.0) go to 10
1305 c write(*,*)'---> ',line
1311 iread=iread+ilen(var)+1
1312 read (line(iread:),*,err=12,end=12) val
1313 c write(*,*)'OK: ',var,' = ',val
1319 #elif (defined CRAY)
1320 integer lennam,lenval,ierror
1322 c getenv using a POSIX call, useful on the T3D
1323 c Sept 1996, comment out error check on advice of H. Pritchard
1326 if(lennam.le.0) stop '--error calling getenv--'
1327 call pxfgetenv(var,lennam,val,lenval,ierror)
1328 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1330 call getenv(var,val)