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)
788 & nterm_kcc(j,i),nterm_kcc_Tb(j,i)
789 C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
790 do k=1,nterm_kcc(j,i)
791 do l=1,nterm_kcc_Tb(j,i)
792 read (itorkcc,*,end=121,err=121)
795 do l=1,nterm_kcc_Tb(j,i)
796 read (itorkcc,*,end=121,err=121)
799 read (itorkcc,*,end=121,err=121)
801 read (itorkcc,*,end=121,err=121)
806 C here will be the apropriate recalibrating for D-aminoacid
807 C read (ithetkcc,*,end=121,err=121) nkcctyp
809 read (ithetkcc,*,end=121,err=121) nbend_kcc_Tb(i)
810 do j=1,nbend_kcc_Tb(i)
811 read (ithetkcc,*,end=121,err=121)
815 C Read of Side-chain backbone correlation parameters
816 C Modified 11 May 2012 by Adasko
819 read (isccor,*,end=119,err=119) nsccortyp
821 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
823 isccortyp(i)=-isccortyp(-i)
825 iscprol=isccortyp(20)
826 c write (iout,*) 'ntortyp',ntortyp
828 cc maxinter is maximum interaction sites
832 read (isccor,*,end=119,err=119)
833 &nterm_sccor(i,j),nlor_sccor(i,j)
839 nterm_sccor(-i,j)=nterm_sccor(i,j)
840 nterm_sccor(-i,-j)=nterm_sccor(i,j)
841 nterm_sccor(i,-j)=nterm_sccor(i,j)
842 do k=1,nterm_sccor(i,j)
843 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
845 if (j.eq.iscprol) then
846 if (i.eq.isccortyp(10)) then
847 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
848 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
850 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
851 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
852 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
853 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
854 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
855 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)
860 if (i.eq.isccortyp(10)) then
861 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
862 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
864 if (j.eq.isccortyp(10)) then
865 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
866 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
868 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
869 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
870 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
871 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
872 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
873 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
877 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
878 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
879 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
880 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
883 do k=1,nlor_sccor(i,j)
884 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
885 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
886 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
887 &(1+vlor3sccor(k,i,j)**2)
889 v0sccor(l,i,j)=v0ijsccor
890 v0sccor(l,-i,j)=v0ijsccor1
891 v0sccor(l,i,-j)=v0ijsccor2
892 v0sccor(l,-i,-j)=v0ijsccor3
898 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
899 c write (iout,*) 'ntortyp',ntortyp
901 cc maxinter is maximum interaction sites
905 read (isccor,*,end=119,err=119)
906 & nterm_sccor(i,j),nlor_sccor(i,j)
910 do k=1,nterm_sccor(i,j)
911 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
913 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
916 do k=1,nlor_sccor(i,j)
917 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
918 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
919 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
920 &(1+vlor3sccor(k,i,j)**2)
922 v0sccor(i,j,iblock)=v0ijsccor
930 write (iout,'(/a/)') 'Torsional constants:'
933 write (iout,*) 'ityp',i,' jtyp',j
934 write (iout,*) 'Fourier constants'
935 do k=1,nterm_sccor(i,j)
936 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
938 write (iout,*) 'Lorenz constants'
939 do k=1,nlor_sccor(i,j)
940 write (iout,'(3(1pe15.5))')
941 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
948 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
949 C interaction energy of the Gly, Ala, and Pro prototypes.
953 write (iout,*) "Coefficients of the cumulants"
955 read (ifourier,*) nloctyp
958 read (ifourier,*,end=115,err=115)
959 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
961 read (ifourier,*,end=115,err=115) (bnew1(ii,1,i),ii=1,3)
962 read (ifourier,*,end=115,err=115) (bnew2(ii,1,i),ii=1,3)
963 read (ifourier,*,end=115,err=115) (bnew1(ii,2,i),ii=1,1)
964 read (ifourier,*,end=115,err=115) (bnew2(ii,2,i),ii=1,1)
965 read (ifourier,*,end=115,err=115) (eenew(ii,i),ii=1,1)
968 write (iout,*) 'Type',i
969 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
977 c B1tilde(1,i) = b(3)
978 c B1tilde(2,i) =-b(5)
979 c B1tilde(1,-i) =-b(3)
980 c B1tilde(2,-i) =b(5)
987 B1tilde(1,i) = b(3,i)
988 B1tilde(2,i) =-b(5,i)
989 C B1tilde(1,-i) =-b(3,i)
990 C B1tilde(2,-i) =b(5,i)
1012 Ctilde(1,1,i)=b(7,i)
1013 Ctilde(1,2,i)=b(9,i)
1014 Ctilde(2,1,i)=-b(9,i)
1015 Ctilde(2,2,i)=b(7,i)
1016 Ctilde(1,1,-i)=b(7,i)
1017 Ctilde(1,2,-i)=-b(9,i)
1018 Ctilde(2,1,-i)=b(9,i)
1019 Ctilde(2,2,-i)=b(7,i)
1021 c Ctilde(1,1,i)=0.0d0
1022 c Ctilde(1,2,i)=0.0d0
1023 c Ctilde(2,1,i)=0.0d0
1024 c Ctilde(2,2,i)=0.0d0
1037 Dtilde(1,1,i)=b(6,i)
1038 Dtilde(1,2,i)=b(8,i)
1039 Dtilde(2,1,i)=-b(8,i)
1040 Dtilde(2,2,i)=b(6,i)
1041 Dtilde(1,1,-i)=b(6,i)
1042 Dtilde(1,2,-i)=-b(8,i)
1043 Dtilde(2,1,-i)=b(8,i)
1044 Dtilde(2,2,-i)=b(6,i)
1046 c Dtilde(1,1,i)=0.0d0
1047 c Dtilde(1,2,i)=0.0d0
1048 c Dtilde(2,1,i)=0.0d0
1049 c Dtilde(2,2,i)=0.0d0
1050 EEold(1,1,i)= b(10,i)+b(11,i)
1051 EEold(2,2,i)=-b(10,i)+b(11,i)
1052 EEold(2,1,i)= b(12,i)-b(13,i)
1053 EEold(1,2,i)= b(12,i)+b(13,i)
1054 EEold(1,1,-i)= b(10,i)+b(11,i)
1055 EEold(2,2,-i)=-b(10,i)+b(11,i)
1056 EEold(2,1,-i)=-b(12,i)+b(13,i)
1057 EEold(1,2,-i)=-b(12,i)-b(13,i)
1058 write(iout,*) "TU DOCHODZE"
1064 c ee(2,1,i)=ee(1,2,i)
1069 write (iout,*) 'Type',i
1071 write(iout,*) B1(1,i),B1(2,i)
1073 write(iout,*) B2(1,i),B2(2,i)
1076 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1080 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1084 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1091 C Read electrostatic-interaction parameters
1095 write (iout,'(/a)') 'Electrostatic interaction constants:'
1096 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1097 & 'IT','JT','APP','BPP','AEL6','AEL3'
1099 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1100 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1101 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1102 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1107 app (i,j)=epp(i,j)*rri*rri
1108 bpp (i,j)=-2.0D0*epp(i,j)*rri
1109 ael6(i,j)=elpp6(i,j)*4.2D0**6
1110 ael3(i,j)=elpp3(i,j)*4.2D0**3
1112 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1113 & ael6(i,j),ael3(i,j)
1118 C Read side-chain interaction parameters.
1120 read (isidep,*,end=117,err=117) ipot,expon
1121 if (ipot.lt.1 .or. ipot.gt.5) then
1122 write (iout,'(2a)') 'Error while reading SC interaction',
1123 & 'potential file - unknown potential type.'
1125 call MPI_Finalize(Ierror)
1131 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1132 & ', exponents are ',expon,2*expon
1133 goto (10,20,30,30,40) ipot
1134 C----------------------- LJ potential ---------------------------------
1135 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1136 & (sigma0(i),i=1,ntyp)
1138 write (iout,'(/a/)') 'Parameters of the LJ 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,a)') 'residue','sigma'
1143 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1146 C----------------------- LJK potential --------------------------------
1147 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1148 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1150 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1151 write (iout,'(a/)') 'The epsilon array:'
1152 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1153 write (iout,'(/a)') 'One-body parameters:'
1154 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1155 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1159 C---------------------- GB or BP potential -----------------------------
1161 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1163 read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
1164 read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
1165 read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
1166 read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
1167 C now we start reading lipid
1169 read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
1170 C print *,"WARNING!!"
1172 C epslip(i,j)=epslip(i,j)+0.05d0
1175 write(iout,*) epslip(1,1),"OK?"
1176 C For the GB potential convert sigma'**2 into chi'
1179 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1183 write (iout,'(/a/)') 'Parameters of the BP potential:'
1184 write (iout,'(a/)') 'The epsilon array:'
1185 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1186 write (iout,'(/a)') 'One-body parameters:'
1187 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1189 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1190 & chip(i),alp(i),i=1,ntyp)
1193 C--------------------- GBV potential -----------------------------------
1194 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1195 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1196 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1198 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1199 write (iout,'(a/)') 'The epsilon array:'
1200 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1201 write (iout,'(/a)') 'One-body parameters:'
1202 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1203 & 's||/s_|_^2',' chip ',' alph '
1204 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1205 & sigii(i),chip(i),alp(i),i=1,ntyp)
1209 C-----------------------------------------------------------------------
1210 C Calculate the "working" parameters of SC interactions.
1214 epslip(i,j)=epslip(j,i)
1219 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1220 sigma(j,i)=sigma(i,j)
1221 rs0(i,j)=dwa16*sigma(i,j)
1225 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1226 & 'Working parameters of the SC interactions:',
1227 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1232 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1241 sigeps=dsign(1.0D0,epsij)
1243 aa_aq(i,j)=epsij*rrij*rrij
1244 bb_aq(i,j)=-sigeps*epsij*rrij
1245 aa_aq(j,i)=aa_aq(i,j)
1246 bb_aq(j,i)=bb_aq(i,j)
1247 epsijlip=epslip(i,j)
1248 sigeps=dsign(1.0D0,epsijlip)
1249 epsijlip=dabs(epsijlip)
1250 aa_lip(i,j)=epsijlip*rrij*rrij
1251 bb_lip(i,j)=-sigeps*epsijlip*rrij
1252 aa_lip(j,i)=aa_lip(i,j)
1253 bb_lip(j,i)=bb_lip(i,j)
1255 sigt1sq=sigma0(i)**2
1256 sigt2sq=sigma0(j)**2
1259 ratsig1=sigt2sq/sigt1sq
1260 ratsig2=1.0D0/ratsig1
1261 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1262 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1263 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1267 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1268 sigmaii(i,j)=rsum_max
1269 sigmaii(j,i)=rsum_max
1271 c sigmaii(i,j)=r0(i,j)
1272 c sigmaii(j,i)=r0(i,j)
1274 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1275 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1276 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1277 augm(i,j)=epsij*r_augm**(2*expon)
1278 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1285 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1286 & restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1287 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1293 C Define the SC-p interaction constants (hard-coded; old style)
1296 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1298 c aad(i,1)=0.3D0*4.0D0**12
1299 C Following line for constants currently implemented
1300 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1301 aad(i,1)=1.5D0*4.0D0**12
1302 c aad(i,1)=0.17D0*5.6D0**12
1304 C "Soft" SC-p repulsion
1306 C Following line for constants currently implemented
1307 c aad(i,1)=0.3D0*4.0D0**6
1308 C "Hard" SC-p repulsion
1309 bad(i,1)=3.0D0*4.0D0**6
1310 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1319 C 8/9/01 Read the SC-p interaction constants from file
1322 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1325 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1326 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1327 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1328 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1332 write (iout,*) "Parameters of SC-p interactions:"
1334 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1335 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1341 C Define the constants of the disulfide bridge
1345 c Old arbitrary potential - commented out.
1350 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1351 c energy surface of diethyl disulfide.
1352 c A. Liwo and U. Kozlowska, 11/24/03
1368 C if(me.eq.king) then
1369 C write (iout,'(/a)') "Disulfide bridge parameters:"
1370 C write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1371 C write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1372 C write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1373 C write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1376 C set the variables used for shielding effect
1377 C write (iout,*) "SHIELD MODE",shield_mode
1378 C if (shield_mode.gt.0) then
1379 C VSolvSphere the volume of solving sphere
1381 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1382 C there will be no distinction between proline peptide group and normal peptide
1383 C group in case of shielding parameters
1384 C VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1385 C VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1386 C write (iout,*) VSolvSphere,VSolvSphere_div
1387 C long axis of side chain
1389 C long_r_sidechain(i)=vbldsc0(1,i)
1390 C short_r_sidechain(i)=sigma0(i)
1392 C lets set the buffor value
1396 111 write (iout,*) "Error reading bending energy parameters."
1398 112 write (iout,*) "Error reading rotamer energy parameters."
1400 113 write (iout,*) "Error reading torsional energy parameters."
1402 114 write (iout,*) "Error reading double torsional energy parameters."
1405 & "Error reading cumulant (multibody energy) parameters."
1407 116 write (iout,*) "Error reading electrostatic energy parameters."
1409 1161 write (iout,*) "Error reading electrostatic energy parameters.Lip"
1411 117 write (iout,*) "Error reading side chain interaction parameters."
1413 118 write (iout,*) "Error reading SCp interaction parameters."
1415 119 write (iout,*) "Error reading SCCOR parameters"
1417 121 write (iout,*) "Error in Czybyshev parameters"
1420 call MPI_Finalize(Ierror)
1427 subroutine getenv_loc(var, val)
1428 character(*) var, val
1431 character(2000) line
1434 open (196,file='env',status='old',readonly,shared)
1436 c write(*,*)'looking for ',var
1437 10 read(196,*,err=11,end=11)line
1438 iread=index(line,var)
1439 c write(*,*)iread,' ',var,' ',line
1440 if (iread.eq.0) go to 10
1441 c write(*,*)'---> ',line
1447 iread=iread+ilen(var)+1
1448 read (line(iread:),*,err=12,end=12) val
1449 c write(*,*)'OK: ',var,' = ',val
1455 #elif (defined CRAY)
1456 integer lennam,lenval,ierror
1458 c getenv using a POSIX call, useful on the T3D
1459 c Sept 1996, comment out error check on advice of H. Pritchard
1462 if(lennam.le.0) stop '--error calling getenv--'
1463 call pxfgetenv(var,lennam,val,lenval,ierror)
1464 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1466 call getenv(var,val)
1468 C set the variables used for shielding effect
1469 C if (shield_mode.gt.0) then
1470 C VSolvSphere the volume of solving sphere
1472 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1473 C there will be no distinction between proline peptide group and normal peptide
1474 C group in case of shielding parameters
1475 C VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1476 C VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1477 C long axis of side chain
1479 C long_r_sidechain(i)=vbldsc0(1,i)
1480 C short_r_sidechain(i)=sigma0(i)
1482 C lets set the buffor value