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 write (iout,*) "iliptranpar",iliptranpar
106 read(iliptranpar,*) pepliptran
108 read(iliptranpar,*) liptranene(i)
113 C Read the parameters of the probability distribution/energy expression
114 C of the virtual-bond valence angles theta
117 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),
118 & (bthet(j,i,1,1),j=1,2)
119 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
120 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
121 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
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)
136 athet(1,i,-1,-1)=athet(1,-i,1,1)
137 athet(2,i,-1,-1)=-athet(2,-i,1,1)
138 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
139 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
140 athet(1,i,-1,1)=athet(1,-i,1,1)
141 athet(2,i,-1,1)=-athet(2,-i,1,1)
142 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
143 bthet(2,i,-1,1)=bthet(2,-i,1,1)
144 athet(1,i,1,-1)=-athet(1,-i,1,1)
145 athet(2,i,1,-1)=athet(2,-i,1,1)
146 bthet(1,i,1,-1)=bthet(1,-i,1,1)
147 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
152 polthet(j,i)=polthet(j,-i)
155 gthet(j,i)=gthet(j,-i)
163 & 'Parameters of the virtual-bond valence angles:'
164 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
165 & ' ATHETA0 ',' A1 ',' A2 ',
168 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
169 & a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
171 write (iout,'(/a/9x,5a/79(1h-))')
172 & 'Parameters of the expression for sigma(theta_c):',
173 & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
174 & ' ALPH3 ',' SIGMA0C '
176 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
177 & (polthet(j,i),j=0,3),sigc0(i)
179 write (iout,'(/a/9x,5a/79(1h-))')
180 & 'Parameters of the second gaussian:',
181 & ' THETA0 ',' SIGMA0 ',' G1 ',
184 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
185 & sig0(i),(gthet(j,i),j=1,3)
189 & 'Parameters of the virtual-bond valence angles:'
190 write (iout,'(/a/9x,5a/79(1h-))')
191 & 'Coefficients of expansion',
192 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
193 & ' b1*10^1 ',' b2*10^1 '
195 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
196 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
197 & (10*bthet(j,i,1,1),j=1,2)
199 write (iout,'(/a/9x,5a/79(1h-))')
200 & 'Parameters of the expression for sigma(theta_c):',
201 & ' alpha0 ',' alph1 ',' alph2 ',
202 & ' alhp3 ',' sigma0c '
204 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
205 & (polthet(j,i),j=0,3),sigc0(i)
207 write (iout,'(/a/9x,5a/79(1h-))')
208 & 'Parameters of the second gaussian:',
209 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
212 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
213 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
219 C Read the parameters of Utheta determined from ab initio surfaces
220 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
222 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,
223 & ntheterm3,nsingle,ndouble
224 write (iout,*) "ithep",ithep
226 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
227 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
229 ithetyp(i)=-ithetyp(-i)
232 do i=-maxthetyp,maxthetyp
233 do j=-maxthetyp,maxthetyp
234 do k=-maxthetyp,maxthetyp
235 aa0thet(i,j,k,iblock)=0.0d0
237 aathet(l,i,j,k,iblock)=0.0d0
241 bbthet(m,l,i,j,k,iblock)=0.0d0
242 ccthet(m,l,i,j,k,iblock)=0.0d0
243 ddthet(m,l,i,j,k,iblock)=0.0d0
244 eethet(m,l,i,j,k,iblock)=0.0d0
250 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
251 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
259 c VAR:iblock means terminally blocking group 1=non-proline 2=proline
261 c VAR:ntethtyp is type of theta potentials type currently 0=glycine
262 c VAR:1=non-glicyne non-proline 2=proline
263 c VAR:negative values for D-aminoacid
265 do j=-nthetyp,nthetyp
266 do k=-nthetyp,nthetyp
267 read (ithep,'(6a)',end=111,err=111) res1
268 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
269 c VAR: aa0thet is variable describing the average value of Foureir
270 c VAR: expansion series
271 c VAR: aathet is foureir expansion in theta/2 angle for full formula
272 c VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
273 Condens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
274 read (ithep,*,end=111,err=111)
275 &(aathet(l,i,j,k,iblock),l=1,ntheterm)
276 read (ithep,*,end=111,err=111)
277 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
278 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
279 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
280 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),
282 read (ithep,*,end=111,err=111)
283 & (((ffthet(llll,lll,ll,i,j,k,iblock),
284 & ffthet(lll,llll,ll,i,j,k,iblock),
285 & ggthet(llll,lll,ll,i,j,k,iblock),
286 & ggthet(lll,llll,ll,i,j,k,iblock),
287 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
292 C For dummy ends assign glycine-type coefficients of theta-only terms; the
293 C coefficients of theta-and-gamma-dependent terms are zero.
294 C IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
295 C RECOMENTDED AFTER VERSION 3.3)
299 c aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
300 c aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
302 c aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
303 c aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
306 c aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
308 c aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
311 C AND COMMENT THE LOOPS BELOW
315 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
316 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
318 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
319 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
322 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
324 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
328 C Substitution for D aminoacids from symmetry.
331 do j=-nthetyp,nthetyp
332 do k=-nthetyp,nthetyp
333 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
335 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
339 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
340 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
341 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
342 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
348 ffthet(llll,lll,ll,i,j,k,iblock)=
349 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
350 ffthet(lll,llll,ll,i,j,k,iblock)=
351 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
352 ggthet(llll,lll,ll,i,j,k,iblock)=
353 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
354 ggthet(lll,llll,ll,i,j,k,iblock)=
355 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
364 C Control printout of the coefficients of virtual-bond-angle potentials
367 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
370 do j=-nthetyp,nthetyp
371 do k=-nthetyp,nthetyp
372 write (iout,'(//4a)')
373 & 'Type ',toronelet(i),toronelet(j),toronelet(k)
374 write (iout,'(//a,10x,a)') " l","a[l]"
375 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
376 write (iout,'(i2,1pe15.5)')
377 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
379 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))')
380 & "b",l,"c",l,"d",l,"e",l
382 write (iout,'(i2,4(1pe15.5))') m,
383 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
384 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
388 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
389 & "f+",l,"f-",l,"g+",l,"g-",l
392 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
393 & ffthet(n,m,l,i,j,k,iblock),
394 & ffthet(m,n,l,i,j,k,iblock),
395 & ggthet(n,m,l,i,j,k,iblock),
396 & ggthet(m,n,l,i,j,k,iblock)
406 write (2,*) "Start reading THETA_PDB",ithep_pdb
409 read (ithep_pdb,*,err=111,end=111)
410 & a0thet(i),(athet(j,i,1,1),j=1,2),
411 & (bthet(j,i,1,1),j=1,2)
412 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
413 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
414 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
418 athet(1,i,1,-1)=athet(1,i,1,1)
419 athet(2,i,1,-1)=athet(2,i,1,1)
420 bthet(1,i,1,-1)=-bthet(1,i,1,1)
421 bthet(2,i,1,-1)=-bthet(2,i,1,1)
422 athet(1,i,-1,1)=-athet(1,i,1,1)
423 athet(2,i,-1,1)=-athet(2,i,1,1)
424 bthet(1,i,-1,1)=bthet(1,i,1,1)
425 bthet(2,i,-1,1)=bthet(2,i,1,1)
429 athet(1,i,-1,-1)=athet(1,-i,1,1)
430 athet(2,i,-1,-1)=-athet(2,-i,1,1)
431 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
432 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
433 athet(1,i,-1,1)=athet(1,-i,1,1)
434 athet(2,i,-1,1)=-athet(2,-i,1,1)
435 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
436 bthet(2,i,-1,1)=bthet(2,-i,1,1)
437 athet(1,i,1,-1)=-athet(1,-i,1,1)
438 athet(2,i,1,-1)=athet(2,-i,1,1)
439 bthet(1,i,1,-1)=bthet(1,-i,1,1)
440 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
445 polthet(j,i)=polthet(j,-i)
448 gthet(j,i)=gthet(j,-i)
451 write (2,*) "End reading THETA_PDB"
457 C Read the parameters of the probability distribution/energy expression
458 C of the side chains.
461 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
465 dsc_inv(i)=1.0D0/dsc(i)
476 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
477 & ((blower(k,l,1),l=1,k),k=1,3)
478 censc(1,1,-i)=censc(1,1,i)
479 censc(2,1,-i)=censc(2,1,i)
480 censc(3,1,-i)=-censc(3,1,i)
482 read (irotam,*,end=112,err=112) bsc(j,i)
483 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
484 & ((blower(k,l,j),l=1,k),k=1,3)
485 censc(1,j,-i)=censc(1,j,i)
486 censc(2,j,-i)=censc(2,j,i)
487 censc(3,j,-i)=-censc(3,j,i)
488 C BSC is amplitude of Gaussian
495 akl=akl+blower(k,m,j)*blower(l,m,j)
499 if (((k.eq.3).and.(l.ne.3))
500 & .or.((l.eq.3).and.(k.ne.3))) then
501 gaussc(k,l,j,-i)=-akl
502 gaussc(l,k,j,-i)=-akl
514 write (iout,'(/a)') 'Parameters of side-chain local geometry'
519 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
520 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
521 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
522 & 'log h',(bsc(j,i),j=1,nlobi)
523 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
524 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
526 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
527 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
530 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
531 write (iout,'(a,f10.4,4(16x,f10.4))')
532 & 'Center ',(bsc(j,i),j=1,nlobi)
533 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
542 C Read scrot parameters for potentials determined from all-atom AM1 calculations
543 C added by Urszula Kozlowska 07/11/2007
546 read (irotam,*,end=112,err=112)
548 read (irotam,*,end=112,err=112)
551 read(irotam,*,end=112,err=112) sc_parmin(j,i)
556 C Read the parameters of the probability distribution/energy expression
557 C of the side chains.
559 write (2,*) "Start reading ROTAM_PDB"
561 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
565 dsc_inv(i)=1.0D0/dsc(i)
576 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
577 & ((blower(k,l,1),l=1,k),k=1,3)
579 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
580 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
581 & ((blower(k,l,j),l=1,k),k=1,3)
588 akl=akl+blower(k,m,j)*blower(l,m,j)
598 write (2,*) "End reading ROTAM_PDB"
604 C Read torsional parameters in old format
606 read (itorp,*,end=113,err=113) ntortyp,nterm_old
607 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
608 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
613 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
619 write (iout,'(/a/)') 'Torsional constants:'
622 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
623 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
629 C Read torsional parameters
631 read (itorp,*,end=113,err=113) ntortyp
632 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
635 itortyp(i)=-itortyp(-i)
637 write (iout,*) 'ntortyp',ntortyp
639 do j=-ntortyp+1,ntortyp-1
640 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
642 nterm(-i,-j,iblock)=nterm(i,j,iblock)
643 nlor(-i,-j,iblock)=nlor(i,j,iblock)
646 do k=1,nterm(i,j,iblock)
647 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
649 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
650 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
651 v0ij=v0ij+si*v1(k,i,j,iblock)
653 c write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
654 c write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
655 c &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
657 do k=1,nlor(i,j,iblock)
658 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
659 & vlor2(k,i,j),vlor3(k,i,j)
660 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
663 v0(-i,-j,iblock)=v0ij
669 write (iout,'(/a/)') 'Torsional constants:'
672 do j=-ntortyp+1,ntortyp-1
673 write (iout,*) 'ityp',i,' jtyp',j
674 write (iout,*) 'Fourier constants'
675 do k=1,nterm(i,j,iblock)
676 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
679 write (iout,*) 'Lorenz constants'
680 do k=1,nlor(i,j,iblock)
681 write (iout,'(3(1pe15.5))')
682 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
690 C 6/23/01 Read parameters for double torsionals
694 do j=-ntortyp+1,ntortyp-1
695 do k=-ntortyp+1,ntortyp-1
696 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
697 c write (iout,*) "OK onelett",
700 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
701 & .or. t3.ne.toronelet(k)) then
702 write (iout,*) "Error in double torsional parameter file",
705 call MPI_Finalize(Ierror)
707 stop "Error in double torsional parameter file"
709 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
710 & ntermd_2(i,j,k,iblock)
711 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
712 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
713 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
714 & ntermd_1(i,j,k,iblock))
715 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
716 & ntermd_1(i,j,k,iblock))
717 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
718 & ntermd_1(i,j,k,iblock))
719 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
720 & ntermd_1(i,j,k,iblock))
721 C Martix of D parameters for one dimesional foureir series
722 do l=1,ntermd_1(i,j,k,iblock)
723 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
724 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
725 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
726 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
727 c write(iout,*) "whcodze" ,
728 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
730 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
731 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
732 & v2s(m,l,i,j,k,iblock),
733 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
734 C Martix of D parameters for two dimesional fourier series
735 do l=1,ntermd_2(i,j,k,iblock)
737 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
738 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
739 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
740 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
749 write (iout,*) 'Constants for double torsionals'
752 do j=-ntortyp+1,ntortyp-1
753 do k=-ntortyp+1,ntortyp-1
754 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
755 & ' nsingle',ntermd_1(i,j,k,iblock),
756 & ' ndouble',ntermd_2(i,j,k,iblock)
758 write (iout,*) 'Single angles:'
759 do l=1,ntermd_1(i,j,k,iblock)
760 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
761 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
762 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
763 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
766 write (iout,*) 'Pairs of angles:'
767 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
768 do l=1,ntermd_2(i,j,k,iblock)
769 write (iout,'(i5,20f10.5)')
770 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
773 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
774 do l=1,ntermd_2(i,j,k,iblock)
775 write (iout,'(i5,20f10.5)')
776 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
777 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
786 C read Czybyshev torsional parameters
787 read (itorkcc,*,end=121,err=121) nkcctyp
788 read (itorkcc,*,end=121,err=121) (itortyp_kcc(i),i=1,ntyp)
790 itortyp_kcc(i)=-itortyp_kcc(-i)
794 C first we read the cos and sin gamma parameters
795 read (itorkcc,*,end=121,err=121)
796 & nterm_kcc(j,i),nterm_kcc_Tb(j,i)
797 C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
798 do k=1,nterm_kcc(j,i)
799 do l=1,nterm_kcc_Tb(j,i)
800 read (itorkcc,*,end=121,err=121) v11_chyb(l,k,j,i)
802 do l=1,nterm_kcc_Tb(j,i)
803 read (itorkcc,*,end=121,err=121) v21_chyb(l,k,j,i)
805 do l=1,nterm_kcc_Tb(j,i)
806 read (itorkcc,*,end=121,err=121) v12_chyb(l,k,j,i)
808 do l=1,nterm_kcc_Tb(j,i)
809 read (itorkcc,*,end=121,err=121) v22_chyb(l,k,j,i)
811 read (itorkcc,*,end=121,err=121) v1_kcc(k,j,i)
812 read (itorkcc,*,end=121,err=121) v2_kcc(k,j,i)
817 c Print valence-torsional parameters
819 & "Parameters of the valence-torsional potentials"
822 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
823 write (iout,'(2a20,a15)') "v_kcc","v1_chyb","v2_chyb"
824 do k=1,nterm_kcc(j,i)
825 write (iout,'(i5,f15.10,i5,2f15.10)')
826 & k,v1_kcc(k,j,i),1,v11_chyb(1,k,j,i),v21_chyb(1,k,j,i)
827 do l=2,nterm_kcc_Tb(j,i)
828 write (iout,'(20x,i5,2f15.10)')
829 & l,v11_chyb(l,k,j,i),v21_chyb(l,k,j,i)
831 write (iout,'(i5,f15.10,i5,2f15.10)')
832 & k,v2_kcc(k,j,i),1,v12_chyb(1,k,j,i),v22_chyb(1,k,j,i)
833 do l=2,nterm_kcc_Tb(j,i)
834 write (iout,'(20x,i5,2f15.10)')
835 & l,v12_chyb(l,k,j,i),v22_chyb(l,k,j,i)
842 C here will be the apropriate recalibrating for D-aminoacid
843 C read (ithetkcc,*,end=121,err=121) nkcctyp
845 read (ithetkcc,*,end=121,err=121) nbend_kcc_Tb(i)
846 do j=1,nbend_kcc_Tb(i)
847 read (ithetkcc,*,end=121,err=121) v1bend_chyb(j,i)
852 & "Parameters of the valence-only potentials"
854 write (iout,'(2a)') "Type ",toronelet(i)
855 do k=1,nbend_kcc_Tb(i)
856 write(iout,'(i5,f15.10)') k,v1bend_chyb(k,i)
860 C Read of Side-chain backbone correlation parameters
861 C Modified 11 May 2012 by Adasko
864 read (isccor,*,end=119,err=119) nsccortyp
866 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
868 isccortyp(i)=-isccortyp(-i)
870 iscprol=isccortyp(20)
871 c write (iout,*) 'ntortyp',ntortyp
873 cc maxinter is maximum interaction sites
877 read (isccor,*,end=119,err=119)
878 &nterm_sccor(i,j),nlor_sccor(i,j)
884 nterm_sccor(-i,j)=nterm_sccor(i,j)
885 nterm_sccor(-i,-j)=nterm_sccor(i,j)
886 nterm_sccor(i,-j)=nterm_sccor(i,j)
887 do k=1,nterm_sccor(i,j)
888 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
890 if (j.eq.iscprol) then
891 if (i.eq.isccortyp(10)) then
892 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
893 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
895 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
896 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
897 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
898 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
899 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
900 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
901 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
902 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
905 if (i.eq.isccortyp(10)) then
906 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
907 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
909 if (j.eq.isccortyp(10)) then
910 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
911 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
913 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
914 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
915 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
916 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
917 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
918 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
922 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
923 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
924 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
925 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
928 do k=1,nlor_sccor(i,j)
929 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
930 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
931 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
932 &(1+vlor3sccor(k,i,j)**2)
934 v0sccor(l,i,j)=v0ijsccor
935 v0sccor(l,-i,j)=v0ijsccor1
936 v0sccor(l,i,-j)=v0ijsccor2
937 v0sccor(l,-i,-j)=v0ijsccor3
943 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
944 c write (iout,*) 'ntortyp',ntortyp
946 cc maxinter is maximum interaction sites
950 read (isccor,*,end=119,err=119)
951 & nterm_sccor(i,j),nlor_sccor(i,j)
955 do k=1,nterm_sccor(i,j)
956 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
958 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
961 do k=1,nlor_sccor(i,j)
962 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
963 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
964 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
965 &(1+vlor3sccor(k,i,j)**2)
967 v0sccor(l,i,j)=v0ijsccor
975 write (iout,'(/a/)') 'Torsional constants:'
979 write (iout,*) 'ityp',i,' jtyp',j
980 write (iout,*) 'Fourier constants'
981 do k=1,nterm_sccor(i,j)
982 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
984 write (iout,*) 'Lorenz constants'
985 do k=1,nlor_sccor(i,j)
986 write (iout,'(3(1pe15.5))')
987 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
995 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
996 C interaction energy of the Gly, Ala, and Pro prototypes.
1000 write (iout,*) "Coefficients of the cumulants"
1002 read (ifourier,*) nloctyp
1004 read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
1005 read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
1006 itype2loc(ntyp1)=nloctyp
1007 iloctyp(nloctyp)=ntyp1
1010 itype2loc(i)=itortyp(i)
1018 itype2loc(-i)=-itype2loc(i)
1021 iloctyp(-i)=-iloctyp(i)
1023 write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1024 write (iout,*) "nloctyp",nloctyp,
1025 & " iloctyp",(iloctyp(i),i=0,nloctyp)
1027 read (ifourier,*,end=115,err=115)
1028 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
1030 read (ifourier,*,end=115,err=115) (bnew1(ii,1,i),ii=1,3)
1031 read (ifourier,*,end=115,err=115) (bnew2(ii,1,i),ii=1,3)
1032 read (ifourier,*,end=115,err=115) (bnew1(ii,2,i),ii=1,1)
1033 read (ifourier,*,end=115,err=115) (bnew2(ii,2,i),ii=1,1)
1034 read (ifourier,*,end=115,err=115) (eenew(ii,i),ii=1,1)
1037 write (iout,*) 'Type',i
1038 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
1046 c B1tilde(1,i) = b(3)
1047 c B1tilde(2,i) =-b(5)
1048 c B1tilde(1,-i) =-b(3)
1049 c B1tilde(2,-i) =b(5)
1050 c b1tilde(1,i)=0.0d0
1051 c b1tilde(2,i)=0.0d0
1056 cc B1tilde(1,i) = b(3,i)
1057 cc B1tilde(2,i) =-b(5,i)
1058 C B1tilde(1,-i) =-b(3,i)
1059 C B1tilde(2,-i) =b(5,i)
1060 cc b1tilde(1,i)=0.0d0
1061 cc b1tilde(2,i)=0.0d0
1081 Ctilde(1,1,i)=b(7,i)
1082 Ctilde(1,2,i)=b(9,i)
1083 Ctilde(2,1,i)=-b(9,i)
1084 Ctilde(2,2,i)=b(7,i)
1085 Ctilde(1,1,-i)=b(7,i)
1086 Ctilde(1,2,-i)=-b(9,i)
1087 Ctilde(2,1,-i)=b(9,i)
1088 Ctilde(2,2,-i)=b(7,i)
1090 c Ctilde(1,1,i)=0.0d0
1091 c Ctilde(1,2,i)=0.0d0
1092 c Ctilde(2,1,i)=0.0d0
1093 c Ctilde(2,2,i)=0.0d0
1106 Dtilde(1,1,i)=b(6,i)
1107 Dtilde(1,2,i)=b(8,i)
1108 Dtilde(2,1,i)=-b(8,i)
1109 Dtilde(2,2,i)=b(6,i)
1110 Dtilde(1,1,-i)=b(6,i)
1111 Dtilde(1,2,-i)=-b(8,i)
1112 Dtilde(2,1,-i)=b(8,i)
1113 Dtilde(2,2,-i)=b(6,i)
1115 c Dtilde(1,1,i)=0.0d0
1116 c Dtilde(1,2,i)=0.0d0
1117 c Dtilde(2,1,i)=0.0d0
1118 c Dtilde(2,2,i)=0.0d0
1119 EEold(1,1,i)= b(10,i)+b(11,i)
1120 EEold(2,2,i)=-b(10,i)+b(11,i)
1121 EEold(2,1,i)= b(12,i)-b(13,i)
1122 EEold(1,2,i)= b(12,i)+b(13,i)
1123 EEold(1,1,-i)= b(10,i)+b(11,i)
1124 EEold(2,2,-i)=-b(10,i)+b(11,i)
1125 EEold(2,1,-i)=-b(12,i)+b(13,i)
1126 EEold(1,2,-i)=-b(12,i)-b(13,i)
1127 write(iout,*) "TU DOCHODZE"
1133 c ee(2,1,i)=ee(1,2,i)
1138 write (iout,*) 'Type',i
1140 write(iout,*) B1(1,i),B1(2,i)
1142 write(iout,*) B2(1,i),B2(2,i)
1145 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1149 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1153 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1160 C Read electrostatic-interaction parameters
1164 write (iout,'(/a)') 'Electrostatic interaction constants:'
1165 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1166 & 'IT','JT','APP','BPP','AEL6','AEL3'
1168 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1169 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1170 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1171 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1176 app (i,j)=epp(i,j)*rri*rri
1177 bpp (i,j)=-2.0D0*epp(i,j)*rri
1178 ael6(i,j)=elpp6(i,j)*4.2D0**6
1179 ael3(i,j)=elpp3(i,j)*4.2D0**3
1181 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1182 & ael6(i,j),ael3(i,j)
1187 C Read side-chain interaction parameters.
1189 read (isidep,*,end=117,err=117) ipot,expon
1190 if (ipot.lt.1 .or. ipot.gt.5) then
1191 write (iout,'(2a)') 'Error while reading SC interaction',
1192 & 'potential file - unknown potential type.'
1194 call MPI_Finalize(Ierror)
1200 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1201 & ', exponents are ',expon,2*expon
1202 goto (10,20,30,30,40) ipot
1203 C----------------------- LJ potential ---------------------------------
1204 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1205 & (sigma0(i),i=1,ntyp)
1207 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1208 write (iout,'(a/)') 'The epsilon array:'
1209 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1210 write (iout,'(/a)') 'One-body parameters:'
1211 write (iout,'(a,4x,a)') 'residue','sigma'
1212 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1215 C----------------------- LJK potential --------------------------------
1216 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1217 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1219 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1220 write (iout,'(a/)') 'The epsilon array:'
1221 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1222 write (iout,'(/a)') 'One-body parameters:'
1223 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1224 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1228 C---------------------- GB or BP potential -----------------------------
1230 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1232 read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
1233 read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
1234 read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
1235 read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
1236 C now we start reading lipid
1238 read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
1239 C print *,"WARNING!!"
1241 C epslip(i,j)=epslip(i,j)+0.05d0
1244 write(iout,*) epslip(1,1),"OK?"
1245 C For the GB potential convert sigma'**2 into chi'
1248 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1252 write (iout,'(/a/)') 'Parameters of the BP potential:'
1253 write (iout,'(a/)') 'The epsilon array:'
1254 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1255 write (iout,'(/a)') 'One-body parameters:'
1256 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1258 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1259 & chip(i),alp(i),i=1,ntyp)
1262 C--------------------- GBV potential -----------------------------------
1263 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1264 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1265 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1267 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1268 write (iout,'(a/)') 'The epsilon array:'
1269 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1270 write (iout,'(/a)') 'One-body parameters:'
1271 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1272 & 's||/s_|_^2',' chip ',' alph '
1273 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1274 & sigii(i),chip(i),alp(i),i=1,ntyp)
1278 C-----------------------------------------------------------------------
1279 C Calculate the "working" parameters of SC interactions.
1283 epslip(i,j)=epslip(j,i)
1288 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1289 sigma(j,i)=sigma(i,j)
1290 rs0(i,j)=dwa16*sigma(i,j)
1294 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1295 & 'Working parameters of the SC interactions:',
1296 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1301 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1310 sigeps=dsign(1.0D0,epsij)
1312 aa_aq(i,j)=epsij*rrij*rrij
1313 bb_aq(i,j)=-sigeps*epsij*rrij
1314 aa_aq(j,i)=aa_aq(i,j)
1315 bb_aq(j,i)=bb_aq(i,j)
1316 epsijlip=epslip(i,j)
1317 sigeps=dsign(1.0D0,epsijlip)
1318 epsijlip=dabs(epsijlip)
1319 aa_lip(i,j)=epsijlip*rrij*rrij
1320 bb_lip(i,j)=-sigeps*epsijlip*rrij
1321 aa_lip(j,i)=aa_lip(i,j)
1322 bb_lip(j,i)=bb_lip(i,j)
1324 sigt1sq=sigma0(i)**2
1325 sigt2sq=sigma0(j)**2
1328 ratsig1=sigt2sq/sigt1sq
1329 ratsig2=1.0D0/ratsig1
1330 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1331 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1332 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1336 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1337 sigmaii(i,j)=rsum_max
1338 sigmaii(j,i)=rsum_max
1340 c sigmaii(i,j)=r0(i,j)
1341 c sigmaii(j,i)=r0(i,j)
1343 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1344 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1345 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1346 augm(i,j)=epsij*r_augm**(2*expon)
1347 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1354 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1355 & restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1356 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1362 C Define the SC-p interaction constants (hard-coded; old style)
1365 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1367 c aad(i,1)=0.3D0*4.0D0**12
1368 C Following line for constants currently implemented
1369 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1370 aad(i,1)=1.5D0*4.0D0**12
1371 c aad(i,1)=0.17D0*5.6D0**12
1373 C "Soft" SC-p repulsion
1375 C Following line for constants currently implemented
1376 c aad(i,1)=0.3D0*4.0D0**6
1377 C "Hard" SC-p repulsion
1378 bad(i,1)=3.0D0*4.0D0**6
1379 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1388 C 8/9/01 Read the SC-p interaction constants from file
1391 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1394 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1395 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1396 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1397 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1401 write (iout,*) "Parameters of SC-p interactions:"
1403 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1404 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1410 C Define the constants of the disulfide bridge
1414 c Old arbitrary potential - commented out.
1419 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1420 c energy surface of diethyl disulfide.
1421 c A. Liwo and U. Kozlowska, 11/24/03
1437 C if(me.eq.king) then
1438 C write (iout,'(/a)') "Disulfide bridge parameters:"
1439 C write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1440 C write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1441 C write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1442 C write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1445 C set the variables used for shielding effect
1446 C write (iout,*) "SHIELD MODE",shield_mode
1447 C if (shield_mode.gt.0) then
1448 C VSolvSphere the volume of solving sphere
1450 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1451 C there will be no distinction between proline peptide group and normal peptide
1452 C group in case of shielding parameters
1453 C VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1454 C VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1455 C write (iout,*) VSolvSphere,VSolvSphere_div
1456 C long axis of side chain
1458 C long_r_sidechain(i)=vbldsc0(1,i)
1459 C short_r_sidechain(i)=sigma0(i)
1461 C lets set the buffor value
1465 111 write (iout,*) "Error reading bending energy parameters."
1467 112 write (iout,*) "Error reading rotamer energy parameters."
1469 113 write (iout,*) "Error reading torsional energy parameters."
1471 114 write (iout,*) "Error reading double torsional energy parameters."
1474 & "Error reading cumulant (multibody energy) parameters."
1476 116 write (iout,*) "Error reading electrostatic energy parameters."
1478 1161 write (iout,*) "Error reading electrostatic energy parameters.Lip"
1480 117 write (iout,*) "Error reading side chain interaction parameters."
1482 118 write (iout,*) "Error reading SCp interaction parameters."
1484 119 write (iout,*) "Error reading SCCOR parameters"
1486 121 write (iout,*) "Error in Czybyshev parameters"
1489 call MPI_Finalize(Ierror)
1496 subroutine getenv_loc(var, val)
1497 character(*) var, val
1500 character(2000) line
1503 open (196,file='env',status='old',readonly,shared)
1505 c write(*,*)'looking for ',var
1506 10 read(196,*,err=11,end=11)line
1507 iread=index(line,var)
1508 c write(*,*)iread,' ',var,' ',line
1509 if (iread.eq.0) go to 10
1510 c write(*,*)'---> ',line
1516 iread=iread+ilen(var)+1
1517 read (line(iread:),*,err=12,end=12) val
1518 c write(*,*)'OK: ',var,' = ',val
1524 #elif (defined CRAY)
1525 integer lennam,lenval,ierror
1527 c getenv using a POSIX call, useful on the T3D
1528 c Sept 1996, comment out error check on advice of H. Pritchard
1531 if(lennam.le.0) stop '--error calling getenv--'
1532 call pxfgetenv(var,lennam,val,lenval,ierror)
1533 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1535 call getenv(var,val)
1537 C set the variables used for shielding effect
1538 C if (shield_mode.gt.0) then
1539 C VSolvSphere the volume of solving sphere
1541 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1542 C there will be no distinction between proline peptide group and normal peptide
1543 C group in case of shielding parameters
1544 C VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1545 C VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1546 C long axis of side chain
1548 C long_r_sidechain(i)=vbldsc0(1,i)
1549 C short_r_sidechain(i)=sigma0(i)
1551 C lets set the buffor value