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)
633 itortyp(ntyp1)=ntortyp
636 itortyp(i)=-itortyp(-i)
638 write (iout,*) 'ntortyp',ntortyp
640 do j=-ntortyp+1,ntortyp-1
641 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
643 nterm(-i,-j,iblock)=nterm(i,j,iblock)
644 nlor(-i,-j,iblock)=nlor(i,j,iblock)
647 do k=1,nterm(i,j,iblock)
648 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
650 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
651 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
652 v0ij=v0ij+si*v1(k,i,j,iblock)
654 c write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
655 c write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
656 c &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
658 do k=1,nlor(i,j,iblock)
659 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
660 & vlor2(k,i,j),vlor3(k,i,j)
661 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
664 v0(-i,-j,iblock)=v0ij
670 write (iout,'(/a/)') 'Torsional constants:'
673 do j=-ntortyp+1,ntortyp-1
674 write (iout,*) 'ityp',i,' jtyp',j
675 write (iout,*) 'Fourier constants'
676 do k=1,nterm(i,j,iblock)
677 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
680 write (iout,*) 'Lorenz constants'
681 do k=1,nlor(i,j,iblock)
682 write (iout,'(3(1pe15.5))')
683 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
691 C 6/23/01 Read parameters for double torsionals
695 do j=-ntortyp+1,ntortyp-1
696 do k=-ntortyp+1,ntortyp-1
697 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
698 c write (iout,*) "OK onelett",
701 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
702 & .or. t3.ne.toronelet(k)) then
703 write (iout,*) "Error in double torsional parameter file",
706 call MPI_Finalize(Ierror)
708 stop "Error in double torsional parameter file"
710 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
711 & ntermd_2(i,j,k,iblock)
712 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
713 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
714 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
715 & ntermd_1(i,j,k,iblock))
716 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
717 & ntermd_1(i,j,k,iblock))
718 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
719 & ntermd_1(i,j,k,iblock))
720 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
721 & ntermd_1(i,j,k,iblock))
722 C Martix of D parameters for one dimesional foureir series
723 do l=1,ntermd_1(i,j,k,iblock)
724 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
725 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
726 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
727 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
728 c write(iout,*) "whcodze" ,
729 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
731 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
732 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
733 & v2s(m,l,i,j,k,iblock),
734 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
735 C Martix of D parameters for two dimesional fourier series
736 do l=1,ntermd_2(i,j,k,iblock)
738 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
739 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
740 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
741 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
750 write (iout,*) 'Constants for double torsionals'
753 do j=-ntortyp+1,ntortyp-1
754 do k=-ntortyp+1,ntortyp-1
755 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
756 & ' nsingle',ntermd_1(i,j,k,iblock),
757 & ' ndouble',ntermd_2(i,j,k,iblock)
759 write (iout,*) 'Single angles:'
760 do l=1,ntermd_1(i,j,k,iblock)
761 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
762 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
763 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
764 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
767 write (iout,*) 'Pairs of angles:'
768 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
769 do l=1,ntermd_2(i,j,k,iblock)
770 write (iout,'(i5,20f10.5)')
771 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
774 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
775 do l=1,ntermd_2(i,j,k,iblock)
776 write (iout,'(i5,20f10.5)')
777 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
778 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
787 C read Czybyshev torsional parameters
788 read (itorkcc,*,end=121,err=121) nkcctyp
789 read (itorkcc,*,end=121,err=121) (itortyp_kcc(i),i=1,ntyp)
791 itortyp_kcc(i)=-itortyp_kcc(-i)
795 C first we read the cos and sin gamma parameters
796 read (itorkcc,*,end=121,err=121)
797 & nterm_kcc(j,i),nterm_kcc_Tb(j,i)
798 C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
799 do k=1,nterm_kcc(j,i)
800 do l=1,nterm_kcc_Tb(j,i)
801 read (itorkcc,*,end=121,err=121) v11_chyb(l,k,j,i)
803 do l=1,nterm_kcc_Tb(j,i)
804 read (itorkcc,*,end=121,err=121) v21_chyb(l,k,j,i)
806 do l=1,nterm_kcc_Tb(j,i)
807 read (itorkcc,*,end=121,err=121) v12_chyb(l,k,j,i)
809 do l=1,nterm_kcc_Tb(j,i)
810 read (itorkcc,*,end=121,err=121) v22_chyb(l,k,j,i)
812 read (itorkcc,*,end=121,err=121) v1_kcc(k,j,i)
813 read (itorkcc,*,end=121,err=121) v2_kcc(k,j,i)
818 c Print valence-torsional parameters
820 & "Parameters of the valence-torsional potentials"
823 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
824 write (iout,'(2a20,a15)') "v_kcc","v1_chyb","v2_chyb"
825 do k=1,nterm_kcc(j,i)
826 write (iout,'(i5,f15.10,i5,2f15.10)')
827 & k,v1_kcc(k,j,i),1,v11_chyb(1,k,j,i),v21_chyb(1,k,j,i)
828 do l=2,nterm_kcc_Tb(j,i)
829 write (iout,'(20x,i5,2f15.10)')
830 & l,v11_chyb(l,k,j,i),v21_chyb(l,k,j,i)
832 write (iout,'(i5,f15.10,i5,2f15.10)')
833 & k,v2_kcc(k,j,i),1,v12_chyb(1,k,j,i),v22_chyb(1,k,j,i)
834 do l=2,nterm_kcc_Tb(j,i)
835 write (iout,'(20x,i5,2f15.10)')
836 & l,v12_chyb(l,k,j,i),v22_chyb(l,k,j,i)
843 C here will be the apropriate recalibrating for D-aminoacid
844 C read (ithetkcc,*,end=121,err=121) nkcctyp
846 read (ithetkcc,*,end=121,err=121) nbend_kcc_Tb(i)
847 do j=1,nbend_kcc_Tb(i)
848 read (ithetkcc,*,end=121,err=121) v1bend_chyb(j,i)
853 & "Parameters of the valence-only potentials"
855 write (iout,'(2a)') "Type ",toronelet(i)
856 do k=1,nbend_kcc_Tb(i)
857 write(iout,'(i5,f15.10)') k,v1bend_chyb(k,i)
861 C Read of Side-chain backbone correlation parameters
862 C Modified 11 May 2012 by Adasko
865 read (isccor,*,end=119,err=119) nsccortyp
867 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
869 isccortyp(i)=-isccortyp(-i)
871 iscprol=isccortyp(20)
872 c write (iout,*) 'ntortyp',ntortyp
874 cc maxinter is maximum interaction sites
878 read (isccor,*,end=119,err=119)
879 &nterm_sccor(i,j),nlor_sccor(i,j)
885 nterm_sccor(-i,j)=nterm_sccor(i,j)
886 nterm_sccor(-i,-j)=nterm_sccor(i,j)
887 nterm_sccor(i,-j)=nterm_sccor(i,j)
888 do k=1,nterm_sccor(i,j)
889 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
891 if (j.eq.iscprol) then
892 if (i.eq.isccortyp(10)) then
893 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
894 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
896 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
897 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
898 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
899 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
900 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
901 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
902 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
903 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
906 if (i.eq.isccortyp(10)) then
907 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
908 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
910 if (j.eq.isccortyp(10)) then
911 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
912 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
914 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
915 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
916 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
917 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
918 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
919 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
923 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
924 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
925 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
926 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
929 do k=1,nlor_sccor(i,j)
930 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
931 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
932 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
933 &(1+vlor3sccor(k,i,j)**2)
935 v0sccor(l,i,j)=v0ijsccor
936 v0sccor(l,-i,j)=v0ijsccor1
937 v0sccor(l,i,-j)=v0ijsccor2
938 v0sccor(l,-i,-j)=v0ijsccor3
944 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
945 c write (iout,*) 'ntortyp',ntortyp
947 cc maxinter is maximum interaction sites
951 read (isccor,*,end=119,err=119)
952 & nterm_sccor(i,j),nlor_sccor(i,j)
956 do k=1,nterm_sccor(i,j)
957 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
959 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
962 do k=1,nlor_sccor(i,j)
963 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
964 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
965 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
966 &(1+vlor3sccor(k,i,j)**2)
968 v0sccor(l,i,j)=v0ijsccor
976 write (iout,'(/a/)') 'Torsional constants:'
980 write (iout,*) 'ityp',i,' jtyp',j
981 write (iout,*) 'Fourier constants'
982 do k=1,nterm_sccor(i,j)
983 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
985 write (iout,*) 'Lorenz constants'
986 do k=1,nlor_sccor(i,j)
987 write (iout,'(3(1pe15.5))')
988 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
996 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
997 C interaction energy of the Gly, Ala, and Pro prototypes.
1001 write (iout,*) "Coefficients of the cumulants"
1003 read (ifourier,*) nloctyp
1005 read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
1006 read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
1007 itype2loc(ntyp1)=nloctyp
1008 iloctyp(nloctyp)=ntyp1
1011 itype2loc(i)=itortyp(i)
1019 itype2loc(-i)=-itype2loc(i)
1022 iloctyp(-i)=-iloctyp(i)
1024 write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1025 write (iout,*) "nloctyp",nloctyp,
1026 & " iloctyp",(iloctyp(i),i=0,nloctyp)
1028 read (ifourier,*,end=115,err=115)
1029 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
1031 read (ifourier,*,end=115,err=115) (bnew1(ii,1,i),ii=1,3)
1032 read (ifourier,*,end=115,err=115) (bnew2(ii,1,i),ii=1,3)
1033 read (ifourier,*,end=115,err=115) (bnew1(ii,2,i),ii=1,1)
1034 read (ifourier,*,end=115,err=115) (bnew2(ii,2,i),ii=1,1)
1035 read (ifourier,*,end=115,err=115) (eenew(ii,i),ii=1,1)
1038 write (iout,*) 'Type',i
1039 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
1051 B1tilde(1,i) = b(3,i)
1052 B1tilde(2,i) =-b(5,i)
1053 B1tilde(1,-i) =-b(3,i)
1054 B1tilde(2,-i) =b(5,i)
1055 c b1tilde(1,i)=0.0d0
1056 c b1tilde(2,i)=0.0d0
1061 cc B1tilde(1,i) = b(3,i)
1062 cc B1tilde(2,i) =-b(5,i)
1063 C B1tilde(1,-i) =-b(3,i)
1064 C B1tilde(2,-i) =b(5,i)
1065 cc b1tilde(1,i)=0.0d0
1066 cc b1tilde(2,i)=0.0d0
1086 Ctilde(1,1,i)=b(7,i)
1087 Ctilde(1,2,i)=b(9,i)
1088 Ctilde(2,1,i)=-b(9,i)
1089 Ctilde(2,2,i)=b(7,i)
1090 Ctilde(1,1,-i)=b(7,i)
1091 Ctilde(1,2,-i)=-b(9,i)
1092 Ctilde(2,1,-i)=b(9,i)
1093 Ctilde(2,2,-i)=b(7,i)
1095 c Ctilde(1,1,i)=0.0d0
1096 c Ctilde(1,2,i)=0.0d0
1097 c Ctilde(2,1,i)=0.0d0
1098 c Ctilde(2,2,i)=0.0d0
1111 Dtilde(1,1,i)=b(6,i)
1112 Dtilde(1,2,i)=b(8,i)
1113 Dtilde(2,1,i)=-b(8,i)
1114 Dtilde(2,2,i)=b(6,i)
1115 Dtilde(1,1,-i)=b(6,i)
1116 Dtilde(1,2,-i)=-b(8,i)
1117 Dtilde(2,1,-i)=b(8,i)
1118 Dtilde(2,2,-i)=b(6,i)
1120 c Dtilde(1,1,i)=0.0d0
1121 c Dtilde(1,2,i)=0.0d0
1122 c Dtilde(2,1,i)=0.0d0
1123 c Dtilde(2,2,i)=0.0d0
1124 EEold(1,1,i)= b(10,i)+b(11,i)
1125 EEold(2,2,i)=-b(10,i)+b(11,i)
1126 EEold(2,1,i)= b(12,i)-b(13,i)
1127 EEold(1,2,i)= b(12,i)+b(13,i)
1128 EEold(1,1,-i)= b(10,i)+b(11,i)
1129 EEold(2,2,-i)=-b(10,i)+b(11,i)
1130 EEold(2,1,-i)=-b(12,i)+b(13,i)
1131 EEold(1,2,-i)=-b(12,i)-b(13,i)
1132 write(iout,*) "TU DOCHODZE"
1138 c ee(2,1,i)=ee(1,2,i)
1143 write (iout,*) 'Type',i
1145 write(iout,*) B1(1,i),B1(2,i)
1147 write(iout,*) B2(1,i),B2(2,i)
1150 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1154 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1158 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1165 C Read electrostatic-interaction parameters
1169 write (iout,'(/a)') 'Electrostatic interaction constants:'
1170 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1171 & 'IT','JT','APP','BPP','AEL6','AEL3'
1173 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1174 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1175 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1176 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1181 app (i,j)=epp(i,j)*rri*rri
1182 bpp (i,j)=-2.0D0*epp(i,j)*rri
1183 ael6(i,j)=elpp6(i,j)*4.2D0**6
1184 ael3(i,j)=elpp3(i,j)*4.2D0**3
1186 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1187 & ael6(i,j),ael3(i,j)
1192 C Read side-chain interaction parameters.
1194 read (isidep,*,end=117,err=117) ipot,expon
1195 if (ipot.lt.1 .or. ipot.gt.5) then
1196 write (iout,'(2a)') 'Error while reading SC interaction',
1197 & 'potential file - unknown potential type.'
1199 call MPI_Finalize(Ierror)
1205 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1206 & ', exponents are ',expon,2*expon
1207 goto (10,20,30,30,40) ipot
1208 C----------------------- LJ potential ---------------------------------
1209 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1210 & (sigma0(i),i=1,ntyp)
1212 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1213 write (iout,'(a/)') 'The epsilon array:'
1214 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1215 write (iout,'(/a)') 'One-body parameters:'
1216 write (iout,'(a,4x,a)') 'residue','sigma'
1217 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1220 C----------------------- LJK potential --------------------------------
1221 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1222 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1224 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1225 write (iout,'(a/)') 'The epsilon array:'
1226 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1227 write (iout,'(/a)') 'One-body parameters:'
1228 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1229 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1233 C---------------------- GB or BP potential -----------------------------
1235 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1237 read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
1238 read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
1239 read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
1240 read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
1241 C now we start reading lipid
1243 read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
1244 write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp)
1246 C print *,"WARNING!!"
1248 C epslip(i,j)=epslip(i,j)+0.05d0
1251 write(iout,*) epslip(1,1),"OK?"
1252 C For the GB potential convert sigma'**2 into chi'
1255 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1259 write (iout,'(/a/)') 'Parameters of the BP potential:'
1260 write (iout,'(a/)') 'The epsilon array:'
1261 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1262 write (iout,'(/a)') 'One-body parameters:'
1263 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1265 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1266 & chip(i),alp(i),i=1,ntyp)
1269 C--------------------- GBV potential -----------------------------------
1270 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1271 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1272 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1274 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1275 write (iout,'(a/)') 'The epsilon array:'
1276 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1277 write (iout,'(/a)') 'One-body parameters:'
1278 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1279 & 's||/s_|_^2',' chip ',' alph '
1280 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1281 & sigii(i),chip(i),alp(i),i=1,ntyp)
1285 C-----------------------------------------------------------------------
1286 C Calculate the "working" parameters of SC interactions.
1290 epslip(i,j)=epslip(j,i)
1295 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1296 sigma(j,i)=sigma(i,j)
1297 rs0(i,j)=dwa16*sigma(i,j)
1301 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1302 & 'Working parameters of the SC interactions:',
1303 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1308 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1317 sigeps=dsign(1.0D0,epsij)
1319 aa_aq(i,j)=epsij*rrij*rrij
1320 bb_aq(i,j)=-sigeps*epsij*rrij
1321 aa_aq(j,i)=aa_aq(i,j)
1322 bb_aq(j,i)=bb_aq(i,j)
1323 epsijlip=epslip(i,j)
1324 sigeps=dsign(1.0D0,epsijlip)
1325 epsijlip=dabs(epsijlip)
1326 aa_lip(i,j)=epsijlip*rrij*rrij
1327 bb_lip(i,j)=-sigeps*epsijlip*rrij
1328 aa_lip(j,i)=aa_lip(i,j)
1329 bb_lip(j,i)=bb_lip(i,j)
1331 sigt1sq=sigma0(i)**2
1332 sigt2sq=sigma0(j)**2
1335 ratsig1=sigt2sq/sigt1sq
1336 ratsig2=1.0D0/ratsig1
1337 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1338 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1339 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1343 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1344 sigmaii(i,j)=rsum_max
1345 sigmaii(j,i)=rsum_max
1347 c sigmaii(i,j)=r0(i,j)
1348 c sigmaii(j,i)=r0(i,j)
1350 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1351 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1352 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1353 augm(i,j)=epsij*r_augm**(2*expon)
1354 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1361 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1362 & restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1363 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1367 write(iout,*) "tube param"
1368 read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep,
1369 & ccavtubpep,dcavtubpep,tubetranenepep
1370 sigmapeptube=sigmapeptube**6
1371 sigeps=dsign(1.0D0,epspeptube)
1372 epspeptube=dabs(epspeptube)
1373 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
1374 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
1375 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
1377 read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i),
1378 & ccavtub(i),dcavtub(i),tubetranene(i)
1379 sigmasctube=sigmasctube**6
1380 sigeps=dsign(1.0D0,epssctube)
1381 epssctube=dabs(epssctube)
1382 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
1383 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
1384 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
1389 C Define the SC-p interaction constants (hard-coded; old style)
1392 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1394 c aad(i,1)=0.3D0*4.0D0**12
1395 C Following line for constants currently implemented
1396 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1397 aad(i,1)=1.5D0*4.0D0**12
1398 c aad(i,1)=0.17D0*5.6D0**12
1400 C "Soft" SC-p repulsion
1402 C Following line for constants currently implemented
1403 c aad(i,1)=0.3D0*4.0D0**6
1404 C "Hard" SC-p repulsion
1405 bad(i,1)=3.0D0*4.0D0**6
1406 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1415 C 8/9/01 Read the SC-p interaction constants from file
1418 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1421 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1422 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1423 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1424 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1428 write (iout,*) "Parameters of SC-p interactions:"
1430 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1431 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1437 C Define the constants of the disulfide bridge
1441 c Old arbitrary potential - commented out.
1446 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1447 c energy surface of diethyl disulfide.
1448 c A. Liwo and U. Kozlowska, 11/24/03
1464 C if(me.eq.king) then
1465 C write (iout,'(/a)') "Disulfide bridge parameters:"
1466 C write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1467 C write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1468 C write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1469 C write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1472 C set the variables used for shielding effect
1473 C write (iout,*) "SHIELD MODE",shield_mode
1474 C if (shield_mode.gt.0) then
1475 C VSolvSphere the volume of solving sphere
1477 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1478 C there will be no distinction between proline peptide group and normal peptide
1479 C group in case of shielding parameters
1480 C VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1481 C VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1482 C write (iout,*) VSolvSphere,VSolvSphere_div
1483 C long axis of side chain
1485 C long_r_sidechain(i)=vbldsc0(1,i)
1486 C short_r_sidechain(i)=sigma0(i)
1488 C lets set the buffor value
1492 111 write (iout,*) "Error reading bending energy parameters."
1494 112 write (iout,*) "Error reading rotamer energy parameters."
1496 113 write (iout,*) "Error reading torsional energy parameters."
1498 114 write (iout,*) "Error reading double torsional energy parameters."
1501 & "Error reading cumulant (multibody energy) parameters."
1503 116 write (iout,*) "Error reading electrostatic energy parameters."
1505 1161 write (iout,*) "Error reading electrostatic energy parameters.Lip"
1507 117 write (iout,*) "Error reading side chain interaction parameters."
1509 118 write (iout,*) "Error reading SCp interaction parameters."
1511 119 write (iout,*) "Error reading SCCOR parameters"
1513 121 write (iout,*) "Error in Czybyshev parameters"
1516 call MPI_Finalize(Ierror)
1523 subroutine getenv_loc(var, val)
1524 character(*) var, val
1527 character(2000) line
1530 open (196,file='env',status='old',readonly,shared)
1532 c write(*,*)'looking for ',var
1533 10 read(196,*,err=11,end=11)line
1534 iread=index(line,var)
1535 c write(*,*)iread,' ',var,' ',line
1536 if (iread.eq.0) go to 10
1537 c write(*,*)'---> ',line
1543 iread=iread+ilen(var)+1
1544 read (line(iread:),*,err=12,end=12) val
1545 c write(*,*)'OK: ',var,' = ',val
1551 #elif (defined CRAY)
1552 integer lennam,lenval,ierror
1554 c getenv using a POSIX call, useful on the T3D
1555 c Sept 1996, comment out error check on advice of H. Pritchard
1558 if(lennam.le.0) stop '--error calling getenv--'
1559 call pxfgetenv(var,lennam,val,lenval,ierror)
1560 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1562 call getenv(var,val)
1564 C set the variables used for shielding effect
1565 C if (shield_mode.gt.0) then
1566 C VSolvSphere the volume of solving sphere
1568 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1569 C there will be no distinction between proline peptide group and normal peptide
1570 C group in case of shielding parameters
1571 C VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1572 C VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1573 C long axis of side chain
1575 C long_r_sidechain(i)=vbldsc0(1,i)
1576 C short_r_sidechain(i)=sigma0(i)
1578 C lets set the buffor value