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 isccortyp(i)=-isccortyp(-i)
950 cc maxinter is maximum interaction sites
951 C NOW READING FOR LL aminoacid chirality
955 read (isccor,*,end=119,err=119)
956 & nterm_sccor(i,j),nlor_sccor(i,j)
957 nterm_sccor(-i,-j)=nterm_sccor(i,j)
961 do k=1,nterm_sccor(i,j)
962 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
964 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
966 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
967 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
970 do k=1,nlor_sccor(i,j)
971 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
972 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
973 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
974 &(1+vlor3sccor(k,i,j)**2)
976 v0sccor(l,i,j)=v0ijsccor
977 v0sccor(l,-i,-j)=v0ijsccor
981 C NOW READING FOR LD aminoacid chirality
985 read (isccor,*,end=119,err=119)
986 & nterm_sccor(i,j),nlor_sccor(i,j)
987 nterm_sccor(-i,-j)=nterm_sccor(i,j)
991 do k=1,nterm_sccor(i,j)
992 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
994 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
996 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
997 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1000 do k=1,nlor_sccor(i,j)
1001 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
1002 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1003 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
1004 &(1+vlor3sccor(k,i,j)**2)
1006 v0sccor(l,i,j)=v0ijsccor
1007 v0sccor(l,-i,-j)=v0ijsccor
1015 do k=1,nterm_sccor(i,j)
1016 v1sccor(k,l,10,j)=v1sccor(k,l,10,-j)
1017 v2sccor(k,l,10,j)=-v2sccor(k,l,10,-j)
1018 v1sccor(k,l,j,10)=v1sccor(k,l,-j,10)
1019 v2sccor(k,l,j,10)=-v2sccor(k,l,-j,10)
1028 write (iout,'(/a/)') 'Torsional constants:'
1032 write (iout,*) 'ityp',i,' jtyp',j
1033 write (iout,*) 'Fourier constants'
1034 do k=1,nterm_sccor(i,j)
1035 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1037 write (iout,*) 'Lorenz constants'
1038 do k=1,nlor_sccor(i,j)
1039 write (iout,'(3(1pe15.5))')
1040 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1048 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1049 C interaction energy of the Gly, Ala, and Pro prototypes.
1053 write (iout,*) "Coefficients of the cumulants"
1055 read (ifourier,*) nloctyp
1057 read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
1058 read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
1059 itype2loc(ntyp1)=nloctyp
1060 iloctyp(nloctyp)=ntyp1
1063 itype2loc(i)=itortyp(i)
1071 itype2loc(-i)=-itype2loc(i)
1074 iloctyp(-i)=-iloctyp(i)
1076 write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1077 write (iout,*) "nloctyp",nloctyp,
1078 & " iloctyp",(iloctyp(i),i=0,nloctyp)
1080 read (ifourier,*,end=115,err=115)
1081 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
1083 read (ifourier,*,end=115,err=115) (bnew1(ii,1,i),ii=1,3)
1084 read (ifourier,*,end=115,err=115) (bnew2(ii,1,i),ii=1,3)
1085 read (ifourier,*,end=115,err=115) (bnew1(ii,2,i),ii=1,1)
1086 read (ifourier,*,end=115,err=115) (bnew2(ii,2,i),ii=1,1)
1087 read (ifourier,*,end=115,err=115) (eenew(ii,i),ii=1,1)
1090 write (iout,*) 'Type',i
1091 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
1103 B1tilde(1,i) = b(3,i)
1104 B1tilde(2,i) =-b(5,i)
1105 B1tilde(1,-i) =-b(3,i)
1106 B1tilde(2,-i) =b(5,i)
1107 c b1tilde(1,i)=0.0d0
1108 c b1tilde(2,i)=0.0d0
1113 cc B1tilde(1,i) = b(3,i)
1114 cc B1tilde(2,i) =-b(5,i)
1115 C B1tilde(1,-i) =-b(3,i)
1116 C B1tilde(2,-i) =b(5,i)
1117 cc b1tilde(1,i)=0.0d0
1118 cc b1tilde(2,i)=0.0d0
1138 Ctilde(1,1,i)=b(7,i)
1139 Ctilde(1,2,i)=b(9,i)
1140 Ctilde(2,1,i)=-b(9,i)
1141 Ctilde(2,2,i)=b(7,i)
1142 Ctilde(1,1,-i)=b(7,i)
1143 Ctilde(1,2,-i)=-b(9,i)
1144 Ctilde(2,1,-i)=b(9,i)
1145 Ctilde(2,2,-i)=b(7,i)
1147 c Ctilde(1,1,i)=0.0d0
1148 c Ctilde(1,2,i)=0.0d0
1149 c Ctilde(2,1,i)=0.0d0
1150 c Ctilde(2,2,i)=0.0d0
1163 Dtilde(1,1,i)=b(6,i)
1164 Dtilde(1,2,i)=b(8,i)
1165 Dtilde(2,1,i)=-b(8,i)
1166 Dtilde(2,2,i)=b(6,i)
1167 Dtilde(1,1,-i)=b(6,i)
1168 Dtilde(1,2,-i)=-b(8,i)
1169 Dtilde(2,1,-i)=b(8,i)
1170 Dtilde(2,2,-i)=b(6,i)
1172 c Dtilde(1,1,i)=0.0d0
1173 c Dtilde(1,2,i)=0.0d0
1174 c Dtilde(2,1,i)=0.0d0
1175 c Dtilde(2,2,i)=0.0d0
1176 EEold(1,1,i)= b(10,i)+b(11,i)
1177 EEold(2,2,i)=-b(10,i)+b(11,i)
1178 EEold(2,1,i)= b(12,i)-b(13,i)
1179 EEold(1,2,i)= b(12,i)+b(13,i)
1180 EEold(1,1,-i)= b(10,i)+b(11,i)
1181 EEold(2,2,-i)=-b(10,i)+b(11,i)
1182 EEold(2,1,-i)=-b(12,i)+b(13,i)
1183 EEold(1,2,-i)=-b(12,i)-b(13,i)
1184 write(iout,*) "TU DOCHODZE"
1190 c ee(2,1,i)=ee(1,2,i)
1195 write (iout,*) 'Type',i
1197 write(iout,*) B1(1,i),B1(2,i)
1199 write(iout,*) B2(1,i),B2(2,i)
1202 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1206 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1210 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1217 C Read electrostatic-interaction parameters
1221 write (iout,'(/a)') 'Electrostatic interaction constants:'
1222 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1223 & 'IT','JT','APP','BPP','AEL6','AEL3'
1225 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1226 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1227 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1228 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1233 app (i,j)=epp(i,j)*rri*rri
1234 bpp (i,j)=-2.0D0*epp(i,j)*rri
1235 ael6(i,j)=elpp6(i,j)*4.2D0**6
1236 ael3(i,j)=elpp3(i,j)*4.2D0**3
1238 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1239 & ael6(i,j),ael3(i,j)
1244 C Read side-chain interaction parameters.
1246 read (isidep,*,end=117,err=117) ipot,expon
1247 if (ipot.lt.1 .or. ipot.gt.5) then
1248 write (iout,'(2a)') 'Error while reading SC interaction',
1249 & 'potential file - unknown potential type.'
1251 call MPI_Finalize(Ierror)
1257 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1258 & ', exponents are ',expon,2*expon
1259 goto (10,20,30,30,40) ipot
1260 C----------------------- LJ potential ---------------------------------
1261 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1262 & (sigma0(i),i=1,ntyp)
1264 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1265 write (iout,'(a/)') 'The epsilon array:'
1266 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1267 write (iout,'(/a)') 'One-body parameters:'
1268 write (iout,'(a,4x,a)') 'residue','sigma'
1269 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1272 C----------------------- LJK potential --------------------------------
1273 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1274 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1276 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1277 write (iout,'(a/)') 'The epsilon array:'
1278 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1279 write (iout,'(/a)') 'One-body parameters:'
1280 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1281 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1285 C---------------------- GB or BP potential -----------------------------
1287 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1289 read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
1290 read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
1291 read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
1292 read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
1293 C now we start reading lipid
1295 read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
1296 write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp)
1298 C print *,"WARNING!!"
1300 C epslip(i,j)=epslip(i,j)+0.05d0
1303 write(iout,*) epslip(1,1),"OK?"
1304 C For the GB potential convert sigma'**2 into chi'
1307 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1311 write (iout,'(/a/)') 'Parameters of the BP potential:'
1312 write (iout,'(a/)') 'The epsilon array:'
1313 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1314 write (iout,'(/a)') 'One-body parameters:'
1315 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1317 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1318 & chip(i),alp(i),i=1,ntyp)
1321 C--------------------- GBV potential -----------------------------------
1322 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1323 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1324 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1326 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1327 write (iout,'(a/)') 'The epsilon array:'
1328 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1329 write (iout,'(/a)') 'One-body parameters:'
1330 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1331 & 's||/s_|_^2',' chip ',' alph '
1332 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1333 & sigii(i),chip(i),alp(i),i=1,ntyp)
1337 C-----------------------------------------------------------------------
1338 C Calculate the "working" parameters of SC interactions.
1342 epslip(i,j)=epslip(j,i)
1347 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1348 sigma(j,i)=sigma(i,j)
1349 rs0(i,j)=dwa16*sigma(i,j)
1353 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1354 & 'Working parameters of the SC interactions:',
1355 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1360 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1369 sigeps=dsign(1.0D0,epsij)
1371 aa_aq(i,j)=epsij*rrij*rrij
1372 bb_aq(i,j)=-sigeps*epsij*rrij
1373 aa_aq(j,i)=aa_aq(i,j)
1374 bb_aq(j,i)=bb_aq(i,j)
1375 epsijlip=epslip(i,j)
1376 sigeps=dsign(1.0D0,epsijlip)
1377 epsijlip=dabs(epsijlip)
1378 aa_lip(i,j)=epsijlip*rrij*rrij
1379 bb_lip(i,j)=-sigeps*epsijlip*rrij
1380 aa_lip(j,i)=aa_lip(i,j)
1381 bb_lip(j,i)=bb_lip(i,j)
1383 sigt1sq=sigma0(i)**2
1384 sigt2sq=sigma0(j)**2
1387 ratsig1=sigt2sq/sigt1sq
1388 ratsig2=1.0D0/ratsig1
1389 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1390 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1391 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1395 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1396 sigmaii(i,j)=rsum_max
1397 sigmaii(j,i)=rsum_max
1399 c sigmaii(i,j)=r0(i,j)
1400 c sigmaii(j,i)=r0(i,j)
1402 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1403 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1404 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1405 augm(i,j)=epsij*r_augm**(2*expon)
1406 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1413 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1414 & restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1415 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1419 write(iout,*) "tube param"
1420 read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep,
1421 & ccavtubpep,dcavtubpep,tubetranenepep
1422 sigmapeptube=sigmapeptube**6
1423 sigeps=dsign(1.0D0,epspeptube)
1424 epspeptube=dabs(epspeptube)
1425 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
1426 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
1427 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
1429 read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i),
1430 & ccavtub(i),dcavtub(i),tubetranene(i)
1431 sigmasctube=sigmasctube**6
1432 sigeps=dsign(1.0D0,epssctube)
1433 epssctube=dabs(epssctube)
1434 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
1435 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
1436 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
1441 C Define the SC-p interaction constants (hard-coded; old style)
1444 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1446 c aad(i,1)=0.3D0*4.0D0**12
1447 C Following line for constants currently implemented
1448 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1449 aad(i,1)=1.5D0*4.0D0**12
1450 c aad(i,1)=0.17D0*5.6D0**12
1452 C "Soft" SC-p repulsion
1454 C Following line for constants currently implemented
1455 c aad(i,1)=0.3D0*4.0D0**6
1456 C "Hard" SC-p repulsion
1457 bad(i,1)=3.0D0*4.0D0**6
1458 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1467 C 8/9/01 Read the SC-p interaction constants from file
1470 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1473 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1474 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1475 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1476 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1480 write (iout,*) "Parameters of SC-p interactions:"
1482 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1483 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1489 C Define the constants of the disulfide bridge
1493 c Old arbitrary potential - commented out.
1498 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1499 c energy surface of diethyl disulfide.
1500 c A. Liwo and U. Kozlowska, 11/24/03
1516 C if(me.eq.king) then
1517 C write (iout,'(/a)') "Disulfide bridge parameters:"
1518 C write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1519 C write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1520 C write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1521 C write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1524 C set the variables used for shielding effect
1525 C write (iout,*) "SHIELD MODE",shield_mode
1526 C if (shield_mode.gt.0) then
1527 C VSolvSphere the volume of solving sphere
1529 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1530 C there will be no distinction between proline peptide group and normal peptide
1531 C group in case of shielding parameters
1532 C VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1533 C VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1534 C write (iout,*) VSolvSphere,VSolvSphere_div
1535 C long axis of side chain
1537 C long_r_sidechain(i)=vbldsc0(1,i)
1538 C short_r_sidechain(i)=sigma0(i)
1540 C lets set the buffor value
1544 111 write (iout,*) "Error reading bending energy parameters."
1546 112 write (iout,*) "Error reading rotamer energy parameters."
1548 113 write (iout,*) "Error reading torsional energy parameters."
1550 114 write (iout,*) "Error reading double torsional energy parameters."
1553 & "Error reading cumulant (multibody energy) parameters."
1555 116 write (iout,*) "Error reading electrostatic energy parameters."
1557 1161 write (iout,*) "Error reading electrostatic energy parameters.Lip"
1559 117 write (iout,*) "Error reading side chain interaction parameters."
1561 118 write (iout,*) "Error reading SCp interaction parameters."
1563 119 write (iout,*) "Error reading SCCOR parameters"
1565 121 write (iout,*) "Error in Czybyshev parameters"
1568 call MPI_Finalize(Ierror)
1575 subroutine getenv_loc(var, val)
1576 character(*) var, val
1579 character(2000) line
1582 open (196,file='env',status='old',readonly,shared)
1584 c write(*,*)'looking for ',var
1585 10 read(196,*,err=11,end=11)line
1586 iread=index(line,var)
1587 c write(*,*)iread,' ',var,' ',line
1588 if (iread.eq.0) go to 10
1589 c write(*,*)'---> ',line
1595 iread=iread+ilen(var)+1
1596 read (line(iread:),*,err=12,end=12) val
1597 c write(*,*)'OK: ',var,' = ',val
1603 #elif (defined CRAY)
1604 integer lennam,lenval,ierror
1606 c getenv using a POSIX call, useful on the T3D
1607 c Sept 1996, comment out error check on advice of H. Pritchard
1610 if(lennam.le.0) stop '--error calling getenv--'
1611 call pxfgetenv(var,lennam,val,lenval,ierror)
1612 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1614 call getenv(var,val)
1616 C set the variables used for shielding effect
1617 C if (shield_mode.gt.0) then
1618 C VSolvSphere the volume of solving sphere
1620 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1621 C there will be no distinction between proline peptide group and normal peptide
1622 C group in case of shielding parameters
1623 C VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1624 C VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1625 C long axis of side chain
1627 C long_r_sidechain(i)=vbldsc0(1,i)
1628 C short_r_sidechain(i)=sigma0(i)
1630 C lets set the buffor value