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'
29 include 'COMMON.CONTROL'
30 include 'COMMON.SHIELD'
32 character*1 onelett(4) /"G","A","P","D"/
33 character*1 toronelet(-2:2) /"p","a","G","A","P"/
35 dimension blower(3,3,maxlob)
37 character*3 lancuch,ucase
39 C For printing parameters after they are read set the following in the UNRES
42 C setenv PRINT_PARM YES
44 C To print parameters in LaTeX format rather than as ASCII tables:
48 call getenv_loc("PRINT_PARM",lancuch)
49 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
50 call getenv_loc("LATEX",lancuch)
51 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
53 dwa16=2.0d0**(1.0d0/6.0d0)
55 C Assign virtual-bond length
60 c Read the virtual-bond parameters, masses, and moments of inertia
61 c and Stokes' radii of the peptide group and side chains
64 read (ibond,*) vbldp0,vbldpdum,akp,mp,ip,pstok
67 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
72 dsc_inv(i)=1.0D0/dsc(i)
76 read (ibond,*) junk,vbldp0,vbldpdum,akp,rjunk,mp,ip,pstok
79 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
80 & j=1,nbondterm(i)),msc(i),isc(i),restok(i)
85 dsc_inv(i)=1.0D0/dsc(i)
90 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
91 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',
93 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
95 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
96 & vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
98 write (iout,'(13x,3f10.5)')
99 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
103 C reading lipid parameters
104 read(iliptranpar,*) pepliptran
106 read(iliptranpar,*) liptranene(i)
111 C Read the parameters of the probability distribution/energy expression
112 C of the virtual-bond valence angles theta
115 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),
116 & (bthet(j,i,1,1),j=1,2)
117 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
118 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
119 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
123 athet(1,i,1,-1)=athet(1,i,1,1)
124 athet(2,i,1,-1)=athet(2,i,1,1)
125 bthet(1,i,1,-1)=-bthet(1,i,1,1)
126 bthet(2,i,1,-1)=-bthet(2,i,1,1)
127 athet(1,i,-1,1)=-athet(1,i,1,1)
128 athet(2,i,-1,1)=-athet(2,i,1,1)
129 bthet(1,i,-1,1)=bthet(1,i,1,1)
130 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)
138 athet(1,i,-1,1)=athet(1,-i,1,1)
139 athet(2,i,-1,1)=-athet(2,-i,1,1)
140 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
141 bthet(2,i,-1,1)=bthet(2,-i,1,1)
142 athet(1,i,1,-1)=-athet(1,-i,1,1)
143 athet(2,i,1,-1)=athet(2,-i,1,1)
144 bthet(1,i,1,-1)=bthet(1,-i,1,1)
145 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
150 polthet(j,i)=polthet(j,-i)
153 gthet(j,i)=gthet(j,-i)
161 & 'Parameters of the virtual-bond valence angles:'
162 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
163 & ' ATHETA0 ',' A1 ',' A2 ',
166 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
167 & a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
169 write (iout,'(/a/9x,5a/79(1h-))')
170 & 'Parameters of the expression for sigma(theta_c):',
171 & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
172 & ' ALPH3 ',' SIGMA0C '
174 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
175 & (polthet(j,i),j=0,3),sigc0(i)
177 write (iout,'(/a/9x,5a/79(1h-))')
178 & 'Parameters of the second gaussian:',
179 & ' THETA0 ',' SIGMA0 ',' G1 ',
182 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
183 & sig0(i),(gthet(j,i),j=1,3)
187 & 'Parameters of the virtual-bond valence angles:'
188 write (iout,'(/a/9x,5a/79(1h-))')
189 & 'Coefficients of expansion',
190 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
191 & ' b1*10^1 ',' b2*10^1 '
193 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
194 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
195 & (10*bthet(j,i,1,1),j=1,2)
197 write (iout,'(/a/9x,5a/79(1h-))')
198 & 'Parameters of the expression for sigma(theta_c):',
199 & ' alpha0 ',' alph1 ',' alph2 ',
200 & ' alhp3 ',' sigma0c '
202 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
203 & (polthet(j,i),j=0,3),sigc0(i)
205 write (iout,'(/a/9x,5a/79(1h-))')
206 & 'Parameters of the second gaussian:',
207 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
210 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
211 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
217 C Read the parameters of Utheta determined from ab initio surfaces
218 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
220 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,
221 & ntheterm3,nsingle,ndouble
222 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
223 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
225 ithetyp(i)=-ithetyp(-i)
228 do i=-maxthetyp,maxthetyp
229 do j=-maxthetyp,maxthetyp
230 do k=-maxthetyp,maxthetyp
231 aa0thet(i,j,k,iblock)=0.0d0
233 aathet(l,i,j,k,iblock)=0.0d0
237 bbthet(m,l,i,j,k,iblock)=0.0d0
238 ccthet(m,l,i,j,k,iblock)=0.0d0
239 ddthet(m,l,i,j,k,iblock)=0.0d0
240 eethet(m,l,i,j,k,iblock)=0.0d0
246 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
247 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
255 c VAR:iblock means terminally blocking group 1=non-proline 2=proline
257 c VAR:ntethtyp is type of theta potentials type currently 0=glycine
258 c VAR:1=non-glicyne non-proline 2=proline
259 c VAR:negative values for D-aminoacid
261 do j=-nthetyp,nthetyp
262 do k=-nthetyp,nthetyp
263 read (ithep,'(6a)',end=111,err=111) res1
264 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
265 c VAR: aa0thet is variable describing the average value of Foureir
266 c VAR: expansion series
267 c VAR: aathet is foureir expansion in theta/2 angle for full formula
268 c VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
269 Condens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
270 read (ithep,*,end=111,err=111)
271 &(aathet(l,i,j,k,iblock),l=1,ntheterm)
272 read (ithep,*,end=111,err=111)
273 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
274 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
275 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
276 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),
278 read (ithep,*,end=111,err=111)
279 & (((ffthet(llll,lll,ll,i,j,k,iblock),
280 & ffthet(lll,llll,ll,i,j,k,iblock),
281 & ggthet(llll,lll,ll,i,j,k,iblock),
282 & ggthet(lll,llll,ll,i,j,k,iblock),
283 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
288 C For dummy ends assign glycine-type coefficients of theta-only terms; the
289 C coefficients of theta-and-gamma-dependent terms are zero.
290 C IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
291 C RECOMENTDED AFTER VERSION 3.3)
295 c aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
296 c aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
298 c aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
299 c aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
302 c aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
304 c aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
307 C AND COMMENT THE LOOPS BELOW
311 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
312 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
314 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
315 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
318 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
320 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
324 C Substitution for D aminoacids from symmetry.
327 do j=-nthetyp,nthetyp
328 do k=-nthetyp,nthetyp
329 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
331 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
335 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
336 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
337 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
338 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
344 ffthet(llll,lll,ll,i,j,k,iblock)=
345 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
346 ffthet(lll,llll,ll,i,j,k,iblock)=
347 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
348 ggthet(llll,lll,ll,i,j,k,iblock)=
349 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
350 ggthet(lll,llll,ll,i,j,k,iblock)=
351 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
360 C Control printout of the coefficients of virtual-bond-angle potentials
363 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
367 write (iout,'(//4a)')
368 & 'Type ',onelett(i),onelett(j),onelett(k)
369 write (iout,'(//a,10x,a)') " l","a[l]"
370 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
371 write (iout,'(i2,1pe15.5)')
372 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
374 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))')
375 & "b",l,"c",l,"d",l,"e",l
377 write (iout,'(i2,4(1pe15.5))') m,
378 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
379 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
383 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
384 & "f+",l,"f-",l,"g+",l,"g-",l
387 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
388 & ffthet(n,m,l,i,j,k,iblock),
389 & ffthet(m,n,l,i,j,k,iblock),
390 & ggthet(n,m,l,i,j,k,iblock),
391 & ggthet(m,n,l,i,j,k,iblock)
400 write (2,*) "Start reading THETA_PDB",ithep_pdb
403 read (ithep_pdb,*,err=111,end=111)
404 & a0thet(i),(athet(j,i,1,1),j=1,2),
405 & (bthet(j,i,1,1),j=1,2)
406 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
407 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
408 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
412 athet(1,i,1,-1)=athet(1,i,1,1)
413 athet(2,i,1,-1)=athet(2,i,1,1)
414 bthet(1,i,1,-1)=-bthet(1,i,1,1)
415 bthet(2,i,1,-1)=-bthet(2,i,1,1)
416 athet(1,i,-1,1)=-athet(1,i,1,1)
417 athet(2,i,-1,1)=-athet(2,i,1,1)
418 bthet(1,i,-1,1)=bthet(1,i,1,1)
419 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)
427 athet(1,i,-1,1)=athet(1,-i,1,1)
428 athet(2,i,-1,1)=-athet(2,-i,1,1)
429 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
430 bthet(2,i,-1,1)=bthet(2,-i,1,1)
431 athet(1,i,1,-1)=-athet(1,-i,1,1)
432 athet(2,i,1,-1)=athet(2,-i,1,1)
433 bthet(1,i,1,-1)=bthet(1,-i,1,1)
434 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
439 polthet(j,i)=polthet(j,-i)
442 gthet(j,i)=gthet(j,-i)
445 write (2,*) "End reading THETA_PDB"
451 C Read the parameters of the probability distribution/energy expression
452 C of the side chains.
455 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
459 dsc_inv(i)=1.0D0/dsc(i)
470 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
471 & ((blower(k,l,1),l=1,k),k=1,3)
472 censc(1,1,-i)=censc(1,1,i)
473 censc(2,1,-i)=censc(2,1,i)
474 censc(3,1,-i)=-censc(3,1,i)
476 read (irotam,*,end=112,err=112) bsc(j,i)
477 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
478 & ((blower(k,l,j),l=1,k),k=1,3)
479 censc(1,j,-i)=censc(1,j,i)
480 censc(2,j,-i)=censc(2,j,i)
481 censc(3,j,-i)=-censc(3,j,i)
482 C BSC is amplitude of Gaussian
489 akl=akl+blower(k,m,j)*blower(l,m,j)
493 if (((k.eq.3).and.(l.ne.3))
494 & .or.((l.eq.3).and.(k.ne.3))) then
495 gaussc(k,l,j,-i)=-akl
496 gaussc(l,k,j,-i)=-akl
508 write (iout,'(/a)') 'Parameters of side-chain local geometry'
513 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
514 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
515 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
516 & 'log h',(bsc(j,i),j=1,nlobi)
517 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
518 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
520 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
521 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
524 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
525 write (iout,'(a,f10.4,4(16x,f10.4))')
526 & 'Center ',(bsc(j,i),j=1,nlobi)
527 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
536 C Read scrot parameters for potentials determined from all-atom AM1 calculations
537 C added by Urszula Kozlowska 07/11/2007
540 read (irotam,*,end=112,err=112)
542 read (irotam,*,end=112,err=112)
545 read(irotam,*,end=112,err=112) sc_parmin(j,i)
550 C Read the parameters of the probability distribution/energy expression
551 C of the side chains.
553 write (2,*) "Start reading ROTAM_PDB"
555 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
559 dsc_inv(i)=1.0D0/dsc(i)
570 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
571 & ((blower(k,l,1),l=1,k),k=1,3)
573 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
574 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
575 & ((blower(k,l,j),l=1,k),k=1,3)
582 akl=akl+blower(k,m,j)*blower(l,m,j)
592 write (2,*) "End reading ROTAM_PDB"
598 C Read torsional parameters in old format
600 read (itorp,*,end=113,err=113) ntortyp,nterm_old
601 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
602 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
607 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
613 write (iout,'(/a/)') 'Torsional constants:'
616 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
617 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
623 C Read torsional parameters
625 read (itorp,*,end=113,err=113) ntortyp
626 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
629 itortyp(i)=-itortyp(-i)
631 write (iout,*) 'ntortyp',ntortyp
633 do j=-ntortyp+1,ntortyp-1
634 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
636 nterm(-i,-j,iblock)=nterm(i,j,iblock)
637 nlor(-i,-j,iblock)=nlor(i,j,iblock)
640 do k=1,nterm(i,j,iblock)
641 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
643 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
644 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
645 v0ij=v0ij+si*v1(k,i,j,iblock)
647 c write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
648 c write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
649 c &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
651 do k=1,nlor(i,j,iblock)
652 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
653 & vlor2(k,i,j),vlor3(k,i,j)
654 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
657 v0(-i,-j,iblock)=v0ij
663 write (iout,'(/a/)') 'Torsional constants:'
666 write (iout,*) 'ityp',i,' jtyp',j
667 write (iout,*) 'Fourier constants'
668 do k=1,nterm(i,j,iblock)
669 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
672 write (iout,*) 'Lorenz constants'
673 do k=1,nlor(i,j,iblock)
674 write (iout,'(3(1pe15.5))')
675 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
682 C 6/23/01 Read parameters for double torsionals
686 do j=-ntortyp+1,ntortyp-1
687 do k=-ntortyp+1,ntortyp-1
688 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
689 c write (iout,*) "OK onelett",
692 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
693 & .or. t3.ne.toronelet(k)) then
694 write (iout,*) "Error in double torsional parameter file",
697 call MPI_Finalize(Ierror)
699 stop "Error in double torsional parameter file"
701 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
702 & ntermd_2(i,j,k,iblock)
703 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
704 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
705 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
706 & ntermd_1(i,j,k,iblock))
707 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
708 & ntermd_1(i,j,k,iblock))
709 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
710 & ntermd_1(i,j,k,iblock))
711 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
712 & ntermd_1(i,j,k,iblock))
713 C Martix of D parameters for one dimesional foureir series
714 do l=1,ntermd_1(i,j,k,iblock)
715 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
716 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
717 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
718 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
719 c write(iout,*) "whcodze" ,
720 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
722 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
723 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
724 & v2s(m,l,i,j,k,iblock),
725 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
726 C Martix of D parameters for two dimesional fourier series
727 do l=1,ntermd_2(i,j,k,iblock)
729 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
730 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
731 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
732 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
741 write (iout,*) 'Constants for double torsionals'
744 do j=-ntortyp+1,ntortyp-1
745 do k=-ntortyp+1,ntortyp-1
746 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
747 & ' nsingle',ntermd_1(i,j,k,iblock),
748 & ' ndouble',ntermd_2(i,j,k,iblock)
750 write (iout,*) 'Single angles:'
751 do l=1,ntermd_1(i,j,k,iblock)
752 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
753 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
754 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
755 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
758 write (iout,*) 'Pairs of angles:'
759 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
760 do l=1,ntermd_2(i,j,k,iblock)
761 write (iout,'(i5,20f10.5)')
762 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
765 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
766 do l=1,ntermd_2(i,j,k,iblock)
767 write (iout,'(i5,20f10.5)')
768 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
769 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
778 C read Czybyshev torsional parameters
779 read (itorkcc,*,end=121,err=121) nkcctyp
780 read (itorkcc,*,end=121,err=121) (itortyp_kcc(i),i=1,ntyp)
782 itortyp_kcc(i)=-itortyp_kcc(-i)
786 C first we read the cos and sin gamma parameters
787 read (itorkcc,*,end=121,err=121) nterm_kcc(j,i)
788 do k=1,nterm_kcc(j,i)
789 read (itorkcc,*,end=121,err=121)
790 & v1_kcc(k,j,i),v2_kcc(k,j,i)
792 C now Czybyshev parameters
793 read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
794 do k=1,nterm_kcc_Tb(j,i)
795 read (itorkcc,*,end=121,err=121)
796 & v1_chyb(k,j,i),v2_chyb(k,j,i)
800 C here will be the apropriate recalibrating for D-aminoacid
803 C Read of Side-chain backbone correlation parameters
804 C Modified 11 May 2012 by Adasko
807 read (isccor,*,end=119,err=119) nsccortyp
809 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
811 isccortyp(i)=-isccortyp(-i)
813 iscprol=isccortyp(20)
814 c write (iout,*) 'ntortyp',ntortyp
816 cc maxinter is maximum interaction sites
820 read (isccor,*,end=119,err=119)
821 &nterm_sccor(i,j),nlor_sccor(i,j)
827 nterm_sccor(-i,j)=nterm_sccor(i,j)
828 nterm_sccor(-i,-j)=nterm_sccor(i,j)
829 nterm_sccor(i,-j)=nterm_sccor(i,j)
830 do k=1,nterm_sccor(i,j)
831 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
833 if (j.eq.iscprol) then
834 if (i.eq.isccortyp(10)) then
835 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
836 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
838 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
839 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
840 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
841 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
842 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
843 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
844 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
845 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
848 if (i.eq.isccortyp(10)) then
849 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
850 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
852 if (j.eq.isccortyp(10)) then
853 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
854 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
856 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
857 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
858 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
859 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
860 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
861 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
865 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
866 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
867 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
868 v0ijsccor3=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(l,i,j)=v0ijsccor
878 v0sccor(l,-i,j)=v0ijsccor1
879 v0sccor(l,i,-j)=v0ijsccor2
880 v0sccor(l,-i,-j)=v0ijsccor3
886 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
887 c write (iout,*) 'ntortyp',ntortyp
889 cc maxinter is maximum interaction sites
893 read (isccor,*,end=119,err=119)
894 & nterm_sccor(i,j),nlor_sccor(i,j)
898 do k=1,nterm_sccor(i,j)
899 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
901 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
904 do k=1,nlor_sccor(i,j)
905 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
906 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
907 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
908 &(1+vlor3sccor(k,i,j)**2)
910 v0sccor(i,j,iblock)=v0ijsccor
918 write (iout,'(/a/)') 'Torsional constants:'
921 write (iout,*) 'ityp',i,' jtyp',j
922 write (iout,*) 'Fourier constants'
923 do k=1,nterm_sccor(i,j)
924 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
926 write (iout,*) 'Lorenz constants'
927 do k=1,nlor_sccor(i,j)
928 write (iout,'(3(1pe15.5))')
929 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
936 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
937 C interaction energy of the Gly, Ala, and Pro prototypes.
941 write (iout,*) "Coefficients of the cumulants"
943 read (ifourier,*) nloctyp
946 read (ifourier,*,end=115,err=115)
947 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
949 read (ifourier,*,end=115,err=115) (bnew1(ii,1,i),ii=1,3)
950 read (ifourier,*,end=115,err=115) (bnew2(ii,1,i),ii=1,3)
951 read (ifourier,*,end=115,err=115) (bnew1(ii,2,i),ii=1,1)
952 read (ifourier,*,end=115,err=115) (bnew2(ii,2,i),ii=1,1)
953 read (ifourier,*,end=115,err=115) (eenew(ii,i),ii=1,1)
956 write (iout,*) 'Type',i
957 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
965 c B1tilde(1,i) = b(3)
966 c B1tilde(2,i) =-b(5)
967 c B1tilde(1,-i) =-b(3)
968 c B1tilde(2,-i) =b(5)
975 B1tilde(1,i) = b(3,i)
976 B1tilde(2,i) =-b(5,i)
977 C B1tilde(1,-i) =-b(3,i)
978 C B1tilde(2,-i) =b(5,i)
1000 Ctilde(1,1,i)=b(7,i)
1001 Ctilde(1,2,i)=b(9,i)
1002 Ctilde(2,1,i)=-b(9,i)
1003 Ctilde(2,2,i)=b(7,i)
1004 Ctilde(1,1,-i)=b(7,i)
1005 Ctilde(1,2,-i)=-b(9,i)
1006 Ctilde(2,1,-i)=b(9,i)
1007 Ctilde(2,2,-i)=b(7,i)
1009 c Ctilde(1,1,i)=0.0d0
1010 c Ctilde(1,2,i)=0.0d0
1011 c Ctilde(2,1,i)=0.0d0
1012 c Ctilde(2,2,i)=0.0d0
1025 Dtilde(1,1,i)=b(6,i)
1026 Dtilde(1,2,i)=b(8,i)
1027 Dtilde(2,1,i)=-b(8,i)
1028 Dtilde(2,2,i)=b(6,i)
1029 Dtilde(1,1,-i)=b(6,i)
1030 Dtilde(1,2,-i)=-b(8,i)
1031 Dtilde(2,1,-i)=b(8,i)
1032 Dtilde(2,2,-i)=b(6,i)
1034 c Dtilde(1,1,i)=0.0d0
1035 c Dtilde(1,2,i)=0.0d0
1036 c Dtilde(2,1,i)=0.0d0
1037 c Dtilde(2,2,i)=0.0d0
1038 EEold(1,1,i)= b(10,i)+b(11,i)
1039 EEold(2,2,i)=-b(10,i)+b(11,i)
1040 EEold(2,1,i)= b(12,i)-b(13,i)
1041 EEold(1,2,i)= b(12,i)+b(13,i)
1042 EEold(1,1,-i)= b(10,i)+b(11,i)
1043 EEold(2,2,-i)=-b(10,i)+b(11,i)
1044 EEold(2,1,-i)=-b(12,i)+b(13,i)
1045 EEold(1,2,-i)=-b(12,i)-b(13,i)
1046 write(iout,*) "TU DOCHODZE"
1052 c ee(2,1,i)=ee(1,2,i)
1057 write (iout,*) 'Type',i
1059 write(iout,*) B1(1,i),B1(2,i)
1061 write(iout,*) B2(1,i),B2(2,i)
1064 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1068 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1072 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1079 C Read electrostatic-interaction parameters
1083 write (iout,'(/a)') 'Electrostatic interaction constants:'
1084 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1085 & 'IT','JT','APP','BPP','AEL6','AEL3'
1087 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1088 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1089 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1090 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1095 app (i,j)=epp(i,j)*rri*rri
1096 bpp (i,j)=-2.0D0*epp(i,j)*rri
1097 ael6(i,j)=elpp6(i,j)*4.2D0**6
1098 ael3(i,j)=elpp3(i,j)*4.2D0**3
1100 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1101 & ael6(i,j),ael3(i,j)
1106 C Read side-chain interaction parameters.
1108 read (isidep,*,end=117,err=117) ipot,expon
1109 if (ipot.lt.1 .or. ipot.gt.5) then
1110 write (iout,'(2a)') 'Error while reading SC interaction',
1111 & 'potential file - unknown potential type.'
1113 call MPI_Finalize(Ierror)
1119 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1120 & ', exponents are ',expon,2*expon
1121 goto (10,20,30,30,40) ipot
1122 C----------------------- LJ potential ---------------------------------
1123 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1124 & (sigma0(i),i=1,ntyp)
1126 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1127 write (iout,'(a/)') 'The epsilon array:'
1128 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1129 write (iout,'(/a)') 'One-body parameters:'
1130 write (iout,'(a,4x,a)') 'residue','sigma'
1131 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1134 C----------------------- LJK potential --------------------------------
1135 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1136 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1138 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1139 write (iout,'(a/)') 'The epsilon array:'
1140 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1141 write (iout,'(/a)') 'One-body parameters:'
1142 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1143 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1147 C---------------------- GB or BP potential -----------------------------
1149 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1151 read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
1152 read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
1153 read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
1154 read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
1155 C now we start reading lipid
1157 read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
1158 C print *,"WARNING!!"
1160 C epslip(i,j)=epslip(i,j)+0.05d0
1163 write(iout,*) epslip(1,1),"OK?"
1164 C For the GB potential convert sigma'**2 into chi'
1167 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1171 write (iout,'(/a/)') 'Parameters of the BP potential:'
1172 write (iout,'(a/)') 'The epsilon array:'
1173 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1174 write (iout,'(/a)') 'One-body parameters:'
1175 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1177 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1178 & chip(i),alp(i),i=1,ntyp)
1181 C--------------------- GBV potential -----------------------------------
1182 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1183 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1184 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1186 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1187 write (iout,'(a/)') 'The epsilon array:'
1188 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1189 write (iout,'(/a)') 'One-body parameters:'
1190 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1191 & 's||/s_|_^2',' chip ',' alph '
1192 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1193 & sigii(i),chip(i),alp(i),i=1,ntyp)
1197 C-----------------------------------------------------------------------
1198 C Calculate the "working" parameters of SC interactions.
1202 epslip(i,j)=epslip(j,i)
1207 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1208 sigma(j,i)=sigma(i,j)
1209 rs0(i,j)=dwa16*sigma(i,j)
1213 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1214 & 'Working parameters of the SC interactions:',
1215 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1220 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1229 sigeps=dsign(1.0D0,epsij)
1231 aa_aq(i,j)=epsij*rrij*rrij
1232 bb_aq(i,j)=-sigeps*epsij*rrij
1233 aa_aq(j,i)=aa_aq(i,j)
1234 bb_aq(j,i)=bb_aq(i,j)
1235 epsijlip=epslip(i,j)
1236 sigeps=dsign(1.0D0,epsijlip)
1237 epsijlip=dabs(epsijlip)
1238 aa_lip(i,j)=epsijlip*rrij*rrij
1239 bb_lip(i,j)=-sigeps*epsijlip*rrij
1240 aa_lip(j,i)=aa_lip(i,j)
1241 bb_lip(j,i)=bb_lip(i,j)
1243 sigt1sq=sigma0(i)**2
1244 sigt2sq=sigma0(j)**2
1247 ratsig1=sigt2sq/sigt1sq
1248 ratsig2=1.0D0/ratsig1
1249 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1250 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1251 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1255 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1256 sigmaii(i,j)=rsum_max
1257 sigmaii(j,i)=rsum_max
1259 c sigmaii(i,j)=r0(i,j)
1260 c sigmaii(j,i)=r0(i,j)
1262 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1263 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1264 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1265 augm(i,j)=epsij*r_augm**(2*expon)
1266 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1273 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1274 & restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1275 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1281 C Define the SC-p interaction constants (hard-coded; old style)
1284 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1286 c aad(i,1)=0.3D0*4.0D0**12
1287 C Following line for constants currently implemented
1288 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1289 aad(i,1)=1.5D0*4.0D0**12
1290 c aad(i,1)=0.17D0*5.6D0**12
1292 C "Soft" SC-p repulsion
1294 C Following line for constants currently implemented
1295 c aad(i,1)=0.3D0*4.0D0**6
1296 C "Hard" SC-p repulsion
1297 bad(i,1)=3.0D0*4.0D0**6
1298 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1307 C 8/9/01 Read the SC-p interaction constants from file
1310 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1313 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1314 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1315 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1316 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1320 write (iout,*) "Parameters of SC-p interactions:"
1322 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1323 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1329 C Define the constants of the disulfide bridge
1333 c Old arbitrary potential - commented out.
1338 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1339 c energy surface of diethyl disulfide.
1340 c A. Liwo and U. Kozlowska, 11/24/03
1356 C if(me.eq.king) then
1357 C write (iout,'(/a)') "Disulfide bridge parameters:"
1358 C write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1359 C write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1360 C write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1361 C write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1364 C set the variables used for shielding effect
1365 C write (iout,*) "SHIELD MODE",shield_mode
1366 C if (shield_mode.gt.0) then
1367 C VSolvSphere the volume of solving sphere
1369 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1370 C there will be no distinction between proline peptide group and normal peptide
1371 C group in case of shielding parameters
1372 C VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1373 C VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1374 C write (iout,*) VSolvSphere,VSolvSphere_div
1375 C long axis of side chain
1377 C long_r_sidechain(i)=vbldsc0(1,i)
1378 C short_r_sidechain(i)=sigma0(i)
1380 C lets set the buffor value
1384 111 write (iout,*) "Error reading bending energy parameters."
1386 112 write (iout,*) "Error reading rotamer energy parameters."
1388 113 write (iout,*) "Error reading torsional energy parameters."
1390 114 write (iout,*) "Error reading double torsional energy parameters."
1393 & "Error reading cumulant (multibody energy) parameters."
1395 116 write (iout,*) "Error reading electrostatic energy parameters."
1397 1161 write (iout,*) "Error reading electrostatic energy parameters.Lip"
1399 117 write (iout,*) "Error reading side chain interaction parameters."
1401 118 write (iout,*) "Error reading SCp interaction parameters."
1403 119 write (iout,*) "Error reading SCCOR parameters"
1405 121 write (iout,*) "Error in Czybyshev parameters"
1408 call MPI_Finalize(Ierror)
1415 subroutine getenv_loc(var, val)
1416 character(*) var, val
1419 character(2000) line
1422 open (196,file='env',status='old',readonly,shared)
1424 c write(*,*)'looking for ',var
1425 10 read(196,*,err=11,end=11)line
1426 iread=index(line,var)
1427 c write(*,*)iread,' ',var,' ',line
1428 if (iread.eq.0) go to 10
1429 c write(*,*)'---> ',line
1435 iread=iread+ilen(var)+1
1436 read (line(iread:),*,err=12,end=12) val
1437 c write(*,*)'OK: ',var,' = ',val
1443 #elif (defined CRAY)
1444 integer lennam,lenval,ierror
1446 c getenv using a POSIX call, useful on the T3D
1447 c Sept 1996, comment out error check on advice of H. Pritchard
1450 if(lennam.le.0) stop '--error calling getenv--'
1451 call pxfgetenv(var,lennam,val,lenval,ierror)
1452 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1454 call getenv(var,val)
1456 C set the variables used for shielding effect
1457 C if (shield_mode.gt.0) then
1458 C VSolvSphere the volume of solving sphere
1460 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1461 C there will be no distinction between proline peptide group and normal peptide
1462 C group in case of shielding parameters
1463 C VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1464 C VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1465 C long axis of side chain
1467 C long_r_sidechain(i)=vbldsc0(1,i)
1468 C short_r_sidechain(i)=sigma0(i)
1470 C lets set the buffor value