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)
1047 c B1tilde(1,i) = b(3)
1048 c B1tilde(2,i) =-b(5)
1049 c B1tilde(1,-i) =-b(3)
1050 c B1tilde(2,-i) =b(5)
1051 c b1tilde(1,i)=0.0d0
1052 c b1tilde(2,i)=0.0d0
1057 cc B1tilde(1,i) = b(3,i)
1058 cc B1tilde(2,i) =-b(5,i)
1059 C B1tilde(1,-i) =-b(3,i)
1060 C B1tilde(2,-i) =b(5,i)
1061 cc b1tilde(1,i)=0.0d0
1062 cc b1tilde(2,i)=0.0d0
1082 Ctilde(1,1,i)=b(7,i)
1083 Ctilde(1,2,i)=b(9,i)
1084 Ctilde(2,1,i)=-b(9,i)
1085 Ctilde(2,2,i)=b(7,i)
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)
1091 c Ctilde(1,1,i)=0.0d0
1092 c Ctilde(1,2,i)=0.0d0
1093 c Ctilde(2,1,i)=0.0d0
1094 c Ctilde(2,2,i)=0.0d0
1107 Dtilde(1,1,i)=b(6,i)
1108 Dtilde(1,2,i)=b(8,i)
1109 Dtilde(2,1,i)=-b(8,i)
1110 Dtilde(2,2,i)=b(6,i)
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)
1116 c Dtilde(1,1,i)=0.0d0
1117 c Dtilde(1,2,i)=0.0d0
1118 c Dtilde(2,1,i)=0.0d0
1119 c Dtilde(2,2,i)=0.0d0
1120 EEold(1,1,i)= b(10,i)+b(11,i)
1121 EEold(2,2,i)=-b(10,i)+b(11,i)
1122 EEold(2,1,i)= b(12,i)-b(13,i)
1123 EEold(1,2,i)= b(12,i)+b(13,i)
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 write(iout,*) "TU DOCHODZE"
1134 c ee(2,1,i)=ee(1,2,i)
1139 write (iout,*) 'Type',i
1141 write(iout,*) B1(1,i),B1(2,i)
1143 write(iout,*) B2(1,i),B2(2,i)
1146 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1150 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1154 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1161 C Read electrostatic-interaction parameters
1165 write (iout,'(/a)') 'Electrostatic interaction constants:'
1166 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1167 & 'IT','JT','APP','BPP','AEL6','AEL3'
1169 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1170 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1171 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1172 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1177 app (i,j)=epp(i,j)*rri*rri
1178 bpp (i,j)=-2.0D0*epp(i,j)*rri
1179 ael6(i,j)=elpp6(i,j)*4.2D0**6
1180 ael3(i,j)=elpp3(i,j)*4.2D0**3
1182 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1183 & ael6(i,j),ael3(i,j)
1188 C Read side-chain interaction parameters.
1190 read (isidep,*,end=117,err=117) ipot,expon
1191 if (ipot.lt.1 .or. ipot.gt.5) then
1192 write (iout,'(2a)') 'Error while reading SC interaction',
1193 & 'potential file - unknown potential type.'
1195 call MPI_Finalize(Ierror)
1201 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1202 & ', exponents are ',expon,2*expon
1203 goto (10,20,30,30,40) ipot
1204 C----------------------- LJ potential ---------------------------------
1205 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1206 & (sigma0(i),i=1,ntyp)
1208 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1209 write (iout,'(a/)') 'The epsilon array:'
1210 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1211 write (iout,'(/a)') 'One-body parameters:'
1212 write (iout,'(a,4x,a)') 'residue','sigma'
1213 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1216 C----------------------- LJK potential --------------------------------
1217 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1218 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1220 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1221 write (iout,'(a/)') 'The epsilon array:'
1222 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1223 write (iout,'(/a)') 'One-body parameters:'
1224 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1225 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1229 C---------------------- GB or BP potential -----------------------------
1231 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1233 read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
1234 read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
1235 read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
1236 read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
1237 C now we start reading lipid
1239 read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
1240 write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp)
1242 C print *,"WARNING!!"
1244 C epslip(i,j)=epslip(i,j)+0.05d0
1247 write(iout,*) epslip(1,1),"OK?"
1248 C For the GB potential convert sigma'**2 into chi'
1251 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1255 write (iout,'(/a/)') 'Parameters of the BP potential:'
1256 write (iout,'(a/)') 'The epsilon array:'
1257 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1258 write (iout,'(/a)') 'One-body parameters:'
1259 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1261 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1262 & chip(i),alp(i),i=1,ntyp)
1265 C--------------------- GBV potential -----------------------------------
1266 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1267 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1268 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1270 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1271 write (iout,'(a/)') 'The epsilon array:'
1272 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1273 write (iout,'(/a)') 'One-body parameters:'
1274 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1275 & 's||/s_|_^2',' chip ',' alph '
1276 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1277 & sigii(i),chip(i),alp(i),i=1,ntyp)
1281 C-----------------------------------------------------------------------
1282 C Calculate the "working" parameters of SC interactions.
1286 epslip(i,j)=epslip(j,i)
1291 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1292 sigma(j,i)=sigma(i,j)
1293 rs0(i,j)=dwa16*sigma(i,j)
1297 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1298 & 'Working parameters of the SC interactions:',
1299 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1304 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1313 sigeps=dsign(1.0D0,epsij)
1315 aa_aq(i,j)=epsij*rrij*rrij
1316 bb_aq(i,j)=-sigeps*epsij*rrij
1317 aa_aq(j,i)=aa_aq(i,j)
1318 bb_aq(j,i)=bb_aq(i,j)
1319 epsijlip=epslip(i,j)
1320 sigeps=dsign(1.0D0,epsijlip)
1321 epsijlip=dabs(epsijlip)
1322 aa_lip(i,j)=epsijlip*rrij*rrij
1323 bb_lip(i,j)=-sigeps*epsijlip*rrij
1324 aa_lip(j,i)=aa_lip(i,j)
1325 bb_lip(j,i)=bb_lip(i,j)
1327 sigt1sq=sigma0(i)**2
1328 sigt2sq=sigma0(j)**2
1331 ratsig1=sigt2sq/sigt1sq
1332 ratsig2=1.0D0/ratsig1
1333 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1334 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1335 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1339 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1340 sigmaii(i,j)=rsum_max
1341 sigmaii(j,i)=rsum_max
1343 c sigmaii(i,j)=r0(i,j)
1344 c sigmaii(j,i)=r0(i,j)
1346 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1347 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1348 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1349 augm(i,j)=epsij*r_augm**(2*expon)
1350 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1357 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1358 & restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1359 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1363 write(iout,*) "tube param"
1364 read(itube,*) epspeptube,sigmapeptube
1365 sigmapeptube=sigmapeptube**6
1366 sigeps=dsign(1.0D0,epspeptube)
1367 epspeptube=dabs(epspeptube)
1368 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
1369 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
1370 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
1372 read(itube,*) epssctube,sigmasctube
1373 sigmasctube=sigmasctube**6
1374 sigeps=dsign(1.0D0,epssctube)
1375 epssctube=dabs(epssctube)
1376 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
1377 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
1378 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
1383 C Define the SC-p interaction constants (hard-coded; old style)
1386 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1388 c aad(i,1)=0.3D0*4.0D0**12
1389 C Following line for constants currently implemented
1390 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1391 aad(i,1)=1.5D0*4.0D0**12
1392 c aad(i,1)=0.17D0*5.6D0**12
1394 C "Soft" SC-p repulsion
1396 C Following line for constants currently implemented
1397 c aad(i,1)=0.3D0*4.0D0**6
1398 C "Hard" SC-p repulsion
1399 bad(i,1)=3.0D0*4.0D0**6
1400 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1409 C 8/9/01 Read the SC-p interaction constants from file
1412 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1415 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1416 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1417 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1418 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1422 write (iout,*) "Parameters of SC-p interactions:"
1424 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1425 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1431 C Define the constants of the disulfide bridge
1435 c Old arbitrary potential - commented out.
1440 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1441 c energy surface of diethyl disulfide.
1442 c A. Liwo and U. Kozlowska, 11/24/03
1458 C if(me.eq.king) then
1459 C write (iout,'(/a)') "Disulfide bridge parameters:"
1460 C write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1461 C write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1462 C write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1463 C write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1466 C set the variables used for shielding effect
1467 C write (iout,*) "SHIELD MODE",shield_mode
1468 C if (shield_mode.gt.0) then
1469 C VSolvSphere the volume of solving sphere
1471 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1472 C there will be no distinction between proline peptide group and normal peptide
1473 C group in case of shielding parameters
1474 C VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1475 C VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1476 C write (iout,*) VSolvSphere,VSolvSphere_div
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
1486 111 write (iout,*) "Error reading bending energy parameters."
1488 112 write (iout,*) "Error reading rotamer energy parameters."
1490 113 write (iout,*) "Error reading torsional energy parameters."
1492 114 write (iout,*) "Error reading double torsional energy parameters."
1495 & "Error reading cumulant (multibody energy) parameters."
1497 116 write (iout,*) "Error reading electrostatic energy parameters."
1499 1161 write (iout,*) "Error reading electrostatic energy parameters.Lip"
1501 117 write (iout,*) "Error reading side chain interaction parameters."
1503 118 write (iout,*) "Error reading SCp interaction parameters."
1505 119 write (iout,*) "Error reading SCCOR parameters"
1507 121 write (iout,*) "Error in Czybyshev parameters"
1510 call MPI_Finalize(Ierror)
1517 subroutine getenv_loc(var, val)
1518 character(*) var, val
1521 character(2000) line
1524 open (196,file='env',status='old',readonly,shared)
1526 c write(*,*)'looking for ',var
1527 10 read(196,*,err=11,end=11)line
1528 iread=index(line,var)
1529 c write(*,*)iread,' ',var,' ',line
1530 if (iread.eq.0) go to 10
1531 c write(*,*)'---> ',line
1537 iread=iread+ilen(var)+1
1538 read (line(iread:),*,err=12,end=12) val
1539 c write(*,*)'OK: ',var,' = ',val
1545 #elif (defined CRAY)
1546 integer lennam,lenval,ierror
1548 c getenv using a POSIX call, useful on the T3D
1549 c Sept 1996, comment out error check on advice of H. Pritchard
1552 if(lennam.le.0) stop '--error calling getenv--'
1553 call pxfgetenv(var,lennam,val,lenval,ierror)
1554 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1556 call getenv(var,val)
1558 C set the variables used for shielding effect
1559 C if (shield_mode.gt.0) then
1560 C VSolvSphere the volume of solving sphere
1562 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1563 C there will be no distinction between proline peptide group and normal peptide
1564 C group in case of shielding parameters
1565 C VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1566 C VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1567 C long axis of side chain
1569 C long_r_sidechain(i)=vbldsc0(1,i)
1570 C short_r_sidechain(i)=sigma0(i)
1572 C lets set the buffor value