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'
30 character*1 onelett(4) /"G","A","P","D"/
31 character*1 toronelet(-2:2) /"p","a","G","A","P"/
33 dimension blower(3,3,maxlob)
35 character*3 lancuch,ucase
37 C For printing parameters after they are read set the following in the UNRES
40 C setenv PRINT_PARM YES
42 C To print parameters in LaTeX format rather than as ASCII tables:
46 call getenv_loc("PRINT_PARM",lancuch)
47 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
48 call getenv_loc("LATEX",lancuch)
49 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
51 dwa16=2.0d0**(1.0d0/6.0d0)
53 C Assign virtual-bond length
58 c Read the virtual-bond parameters, masses, and moments of inertia
59 c and Stokes' radii of the peptide group and side chains
62 read (ibond,*) vbldp0,akp,mp,ip,pstok
65 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
70 dsc_inv(i)=1.0D0/dsc(i)
74 read (ibond,*) junk,vbldp0,akp,rjunk,mp,ip,pstok
76 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
77 & j=1,nbondterm(i)),msc(i),isc(i),restok(i)
82 dsc_inv(i)=1.0D0/dsc(i)
87 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
88 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',
90 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
92 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
93 & vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
95 write (iout,'(13x,3f10.5)')
96 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
102 C Read the parameters of the probability distribution/energy expression
103 C of the virtual-bond valence angles theta
106 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),
107 & (bthet(j,i,1,1),j=1,2)
108 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
109 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
110 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
114 athet(1,i,1,-1)=athet(1,i,1,1)
115 athet(2,i,1,-1)=athet(2,i,1,1)
116 bthet(1,i,1,-1)=-bthet(1,i,1,1)
117 bthet(2,i,1,-1)=-bthet(2,i,1,1)
118 athet(1,i,-1,1)=-athet(1,i,1,1)
119 athet(2,i,-1,1)=-athet(2,i,1,1)
120 bthet(1,i,-1,1)=bthet(1,i,1,1)
121 bthet(2,i,-1,1)=bthet(2,i,1,1)
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)
133 athet(1,i,1,-1)=-athet(1,-i,1,1)
134 athet(2,i,1,-1)=athet(2,-i,1,1)
135 bthet(1,i,1,-1)=bthet(1,-i,1,1)
136 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
141 polthet(j,i)=polthet(j,-i)
144 gthet(j,i)=gthet(j,-i)
152 & 'Parameters of the virtual-bond valence angles:'
153 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
154 & ' ATHETA0 ',' A1 ',' A2 ',
157 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
158 & a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
160 write (iout,'(/a/9x,5a/79(1h-))')
161 & 'Parameters of the expression for sigma(theta_c):',
162 & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
163 & ' ALPH3 ',' SIGMA0C '
165 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
166 & (polthet(j,i),j=0,3),sigc0(i)
168 write (iout,'(/a/9x,5a/79(1h-))')
169 & 'Parameters of the second gaussian:',
170 & ' THETA0 ',' SIGMA0 ',' G1 ',
173 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
174 & sig0(i),(gthet(j,i),j=1,3)
178 & 'Parameters of the virtual-bond valence angles:'
179 write (iout,'(/a/9x,5a/79(1h-))')
180 & 'Coefficients of expansion',
181 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
182 & ' b1*10^1 ',' b2*10^1 '
184 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
185 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
186 & (10*bthet(j,i,1,1),j=1,2)
188 write (iout,'(/a/9x,5a/79(1h-))')
189 & 'Parameters of the expression for sigma(theta_c):',
190 & ' alpha0 ',' alph1 ',' alph2 ',
191 & ' alhp3 ',' sigma0c '
193 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
194 & (polthet(j,i),j=0,3),sigc0(i)
196 write (iout,'(/a/9x,5a/79(1h-))')
197 & 'Parameters of the second gaussian:',
198 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
201 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
202 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
208 C Read the parameters of Utheta determined from ab initio surfaces
209 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
211 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,
212 & ntheterm3,nsingle,ndouble
213 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
214 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
216 ithetyp(i)=-ithetyp(-i)
219 do i=-maxthetyp,maxthetyp
220 do j=-maxthetyp,maxthetyp
221 do k=-maxthetyp,maxthetyp
222 aa0thet(i,j,k,iblock)=0.0d0
224 aathet(l,i,j,k,iblock)=0.0d0
228 bbthet(m,l,i,j,k,iblock)=0.0d0
229 ccthet(m,l,i,j,k,iblock)=0.0d0
230 ddthet(m,l,i,j,k,iblock)=0.0d0
231 eethet(m,l,i,j,k,iblock)=0.0d0
237 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
238 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
246 c VAR:iblock means terminally blocking group 1=non-proline 2=proline
248 c VAR:ntethtyp is type of theta potentials type currently 0=glycine
249 c VAR:1=non-glicyne non-proline 2=proline
250 c VAR:negative values for D-aminoacid
252 do j=-nthetyp,nthetyp
253 do k=-nthetyp,nthetyp
254 read (ithep,'(6a)',end=111,err=111) res1
255 c VAR: aa0thet is variable describing the average value of Foureir
256 c VAR: expansion series
257 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
258 c VAR: aathet is foureir expansion in theta/2 angle for full formula
259 c VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
260 Condens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
261 read (ithep,*,end=111,err=111)
262 &(aathet(l,i,j,k,iblock),l=1,ntheterm)
263 read (ithep,*,end=111,err=111)
264 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
265 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
266 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
267 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),
269 read (ithep,*,end=111,err=111)
270 & (((ffthet(llll,lll,ll,i,j,k,iblock),
271 & ffthet(lll,llll,ll,i,j,k,iblock),
272 & ggthet(llll,lll,ll,i,j,k,iblock),
273 & ggthet(lll,llll,ll,i,j,k,iblock),
274 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
281 C For dummy ends assign glycine-type coefficients of theta-only terms; the
282 C coefficients of theta-and-gamma-dependent terms are zero.
283 C IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
284 C RECOMENTDED AFTER VERSION 3.3)
288 c aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
289 c aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
291 c aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
292 c aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
295 c aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
297 c aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
300 C AND COMMENT THE LOOPS BELOW
304 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
305 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
307 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
308 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
311 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
313 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
317 C Substitution for D aminoacids from symmetry.
320 do j=-nthetyp,nthetyp
321 do k=-nthetyp,nthetyp
322 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
324 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
328 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
329 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
330 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
331 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
337 ffthet(llll,lll,ll,i,j,k,iblock)=
338 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
339 ffthet(lll,llll,ll,i,j,k,iblock)=
340 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
341 ggthet(llll,lll,ll,i,j,k,iblock)=
342 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
343 ggthet(lll,llll,ll,i,j,k,iblock)=
344 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
353 C Control printout of the coefficients of virtual-bond-angle potentials
356 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
360 write (iout,'(//4a)')
361 & 'Type ',onelett(i),onelett(j),onelett(k)
362 write (iout,'(//a,10x,a)') " l","a[l]"
363 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
364 write (iout,'(i2,1pe15.5)')
365 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
367 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))')
368 & "b",l,"c",l,"d",l,"e",l
370 write (iout,'(i2,4(1pe15.5))') m,
371 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
372 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
376 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
377 & "f+",l,"f-",l,"g+",l,"g-",l
380 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
381 & ffthet(n,m,l,i,j,k,iblock),
382 & ffthet(m,n,l,i,j,k,iblock),
383 & ggthet(n,m,l,i,j,k,iblock),
384 & ggthet(m,n,l,i,j,k,iblock)
393 write (2,*) "Start reading THETA_PDB"
395 read (ithep_pdb,*,err=111,end=111)
396 & a0thet(i),(athet(j,i,1,1),j=1,2),
397 & (bthet(j,i,1,1),j=1,2)
398 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
399 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
400 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
404 athet(1,i,1,-1)=athet(1,i,1,1)
405 athet(2,i,1,-1)=athet(2,i,1,1)
406 bthet(1,i,1,-1)=-bthet(1,i,1,1)
407 bthet(2,i,1,-1)=-bthet(2,i,1,1)
408 athet(1,i,-1,1)=-athet(1,i,1,1)
409 athet(2,i,-1,1)=-athet(2,i,1,1)
410 bthet(1,i,-1,1)=bthet(1,i,1,1)
411 bthet(2,i,-1,1)=bthet(2,i,1,1)
415 athet(1,i,-1,-1)=athet(1,-i,1,1)
416 athet(2,i,-1,-1)=-athet(2,-i,1,1)
417 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
418 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
419 athet(1,i,-1,1)=athet(1,-i,1,1)
420 athet(2,i,-1,1)=-athet(2,-i,1,1)
421 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
422 bthet(2,i,-1,1)=bthet(2,-i,1,1)
423 athet(1,i,1,-1)=-athet(1,-i,1,1)
424 athet(2,i,1,-1)=athet(2,-i,1,1)
425 bthet(1,i,1,-1)=bthet(1,-i,1,1)
426 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
431 polthet(j,i)=polthet(j,-i)
434 gthet(j,i)=gthet(j,-i)
437 write (2,*) "End reading THETA_PDB"
443 C Read the parameters of the probability distribution/energy expression
444 C of the side chains.
447 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
451 dsc_inv(i)=1.0D0/dsc(i)
462 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
463 & ((blower(k,l,1),l=1,k),k=1,3)
465 read (irotam,*,end=112,err=112) bsc(j,i)
466 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
467 & ((blower(k,l,j),l=1,k),k=1,3)
474 akl=akl+blower(k,m,j)*blower(l,m,j)
485 write (iout,'(/a)') 'Parameters of side-chain local geometry'
490 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
491 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
492 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
493 & 'log h',(bsc(j,i),j=1,nlobi)
494 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
495 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
497 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
498 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
501 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
502 write (iout,'(a,f10.4,4(16x,f10.4))')
503 & 'Center ',(bsc(j,i),j=1,nlobi)
504 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
513 C Read scrot parameters for potentials determined from all-atom AM1 calculations
514 C added by Urszula Kozlowska 07/11/2007
517 read (irotam,*,end=112,err=112)
519 read (irotam,*,end=112,err=112)
522 read(irotam,*,end=112,err=112) sc_parmin(j,i)
527 C Read the parameters of the probability distribution/energy expression
528 C of the side chains.
530 write (2,*) "Start reading ROTAM_PDB"
533 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
537 dsc_inv(i)=1.0D0/dsc(i)
548 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
549 & ((blower(k,l,1),l=1,k),k=1,3)
551 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
552 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
553 & ((blower(k,l,j),l=1,k),k=1,3)
560 akl=akl+blower(k,m,j)*blower(l,m,j)
570 write (2,*) "Ending reading ROTAM_PDB"
576 C Read torsional parameters in old format
578 read (itorp,*,end=113,err=113) ntortyp,nterm_old
579 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
580 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
585 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
591 write (iout,'(/a/)') 'Torsional constants:'
594 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
595 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
601 C Read torsional parameters
603 read (itorp,*,end=113,err=113) ntortyp
604 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
605 c write (iout,*) 'ntortyp',ntortyp
608 read (itorp,*,end=113,err=113) nterm(i,j),nlor(i,j)
612 read (itorp,*,end=113,err=113) kk,v1(k,i,j),v2(k,i,j)
613 v0ij=v0ij+si*v1(k,i,j)
617 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
618 & vlor2(k,i,j),vlor3(k,i,j)
619 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
626 write (iout,'(/a/)') 'Torsional constants:'
629 write (iout,*) 'ityp',i,' jtyp',j
630 write (iout,*) 'Fourier constants'
632 write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
634 write (iout,*) 'Lorenz constants'
636 write (iout,'(3(1pe15.5))')
637 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
643 C 6/23/01 Read parameters for double torsionals
648 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
649 c write (iout,*) "OK onelett",
652 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
653 & .or. t3.ne.toronelet(k)) then
654 write (iout,*) "Error in double torsional parameter file",
657 call MPI_Finalize(Ierror)
659 stop "Error in double torsional parameter file"
661 read (itordp,*,end=114,err=114) ntermd_1(i,j,k),
663 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k),l=1,
665 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k),l=1,
667 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k),l=1,
669 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k),l=1,
671 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k),
672 & v2c(m,l,i,j,k),v2s(l,m,i,j,k),v2s(m,l,i,j,k),
673 & m=1,l-1),l=1,ntermd_2(i,j,k))
679 write (iout,*) 'Constants for double torsionals'
682 do j=-ntortyp+1,ntortyp-1
683 do k=-ntortyp+1,ntortyp-1
684 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
685 & ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k)
687 write (iout,*) 'Single angles:'
688 do l=1,ntermd_1(i,j,k)
689 write (iout,'(i5,2f10.5,5x,2f10.5)') l,
690 & v1c(1,l,i,j,k),v1s(1,l,i,j,k),
691 & v1c(2,l,i,j,k),v1s(2,l,i,j,k)
694 write (iout,*) 'Pairs of angles:'
695 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
696 do l=1,ntermd_2(i,j,k)
697 write (iout,'(i5,20f10.5)')
698 & l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k))
701 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
702 do l=1,ntermd_2(i,j,k)
703 write (iout,'(i5,20f10.5)')
704 & l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k))
713 C Read of Side-chain backbone correlation parameters
714 C Modified 11 May 2012 by Adasko
717 read (isccor,*,end=1113,err=1113) nsccortyp
719 read (isccor,*,end=1113,err=1113) (isccortyp(i),i=1,ntyp)
721 isccortyp(i)=-isccortyp(-i)
723 iscprol=isccortyp(20)
724 c write (iout,*) 'ntortyp',ntortyp
726 cc maxinter is maximum interaction sites
730 read (isccor,*,end=1113,err=1113) nterm_sccor(i,j),
737 nterm_sccor(-i,j)=nterm_sccor(i,j)
738 nterm_sccor(-i,-j)=nterm_sccor(i,j)
739 nterm_sccor(i,-j)=nterm_sccor(i,j)
740 do k=1,nterm_sccor(i,j)
741 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
743 if (j.eq.iscprol) then
744 if (i.eq.isccortyp(10)) then
745 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
746 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
748 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
749 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
750 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
751 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
752 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
753 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
754 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
755 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
758 if (i.eq.isccortyp(10)) then
759 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
760 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
762 if (j.eq.isccortyp(10)) then
763 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
764 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
766 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
767 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
768 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
769 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
770 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
771 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
775 read (isccor,*,end=1113,err=1113) kk,v1sccor(k,l,i,j)
777 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
778 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
779 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
780 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
783 do k=1,nlor_sccor(i,j)
784 read (isccor,*,end=1113,err=1113) kk,vlor1sccor(k,i,j),
785 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
786 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
787 &(1+vlor3sccor(k,i,j)**2)
789 v0sccor(l,i,j)=v0ijsccor
790 v0sccor(l,-i,j)=v0ijsccor1
791 v0sccor(l,i,-j)=v0ijsccor2
792 v0sccor(l,-i,-j)=v0ijsccor3
798 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
799 c write (iout,*) 'ntortyp',ntortyp
801 cc maxinter is maximum interaction sites
805 read (isccor,*,end=119,err=119)
806 & nterm_sccor(i,j),nlor_sccor(i,j)
810 do k=1,nterm_sccor(i,j)
811 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
813 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
816 do k=1,nlor_sccor(i,j)
817 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
818 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
819 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
820 &(1+vlor3sccor(k,i,j)**2)
822 v0sccor(i,j,iblock)=v0ijsccor
830 write (iout,'(/a/)') 'Torsional constants:'
833 write (iout,*) 'ityp',i,' jtyp',j
834 write (iout,*) 'Fourier constants'
835 do k=1,nterm_sccor(i,j)
836 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
838 write (iout,*) 'Lorenz constants'
839 do k=1,nlor_sccor(i,j)
840 write (iout,'(3(1pe15.5))')
841 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
848 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
849 C interaction energy of the Gly, Ala, and Pro prototypes.
853 write (iout,*) "Coefficients of the cumulants"
855 read (ifourier,*) nloctyp
857 read (ifourier,*,end=115,err=115)
858 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
860 write (iout,*) 'Type',i
861 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
889 c Ctilde(1,1,i)=0.0d0
890 c Ctilde(1,2,i)=0.0d0
891 c Ctilde(2,1,i)=0.0d0
892 c Ctilde(2,2,i)=0.0d0
905 c Dtilde(1,1,i)=0.0d0
906 c Dtilde(1,2,i)=0.0d0
907 c Dtilde(2,1,i)=0.0d0
908 c Dtilde(2,2,i)=0.0d0
909 EE(1,1,i)= b(10)+b(11)
910 EE(2,2,i)=-b(10)+b(11)
911 EE(2,1,i)= b(12)-b(13)
912 EE(1,2,i)= b(12)+b(13)
917 c ee(2,1,i)=ee(1,2,i)
921 write (iout,*) 'Type',i
923 write(iout,*) B1(1,i),B1(2,i)
925 write(iout,*) B2(1,i),B2(2,i)
928 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
932 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
936 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
941 C Read electrostatic-interaction parameters
945 write (iout,'(/a)') 'Electrostatic interaction constants:'
946 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
947 & 'IT','JT','APP','BPP','AEL6','AEL3'
949 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
950 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
951 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
952 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
957 app (i,j)=epp(i,j)*rri*rri
958 bpp (i,j)=-2.0D0*epp(i,j)*rri
959 ael6(i,j)=elpp6(i,j)*4.2D0**6
960 ael3(i,j)=elpp3(i,j)*4.2D0**3
961 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
962 & ael6(i,j),ael3(i,j)
966 C Read side-chain interaction parameters.
968 read (isidep,*,end=117,err=117) ipot,expon
969 if (ipot.lt.1 .or. ipot.gt.5) then
970 write (iout,'(2a)') 'Error while reading SC interaction',
971 & 'potential file - unknown potential type.'
973 call MPI_Finalize(Ierror)
979 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
980 & ', exponents are ',expon,2*expon
981 goto (10,20,30,30,40) ipot
982 C----------------------- LJ potential ---------------------------------
983 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
984 & (sigma0(i),i=1,ntyp)
986 write (iout,'(/a/)') 'Parameters of the LJ potential:'
987 write (iout,'(a/)') 'The epsilon array:'
988 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
989 write (iout,'(/a)') 'One-body parameters:'
990 write (iout,'(a,4x,a)') 'residue','sigma'
991 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
994 C----------------------- LJK potential --------------------------------
995 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
996 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
998 write (iout,'(/a/)') 'Parameters of the LJK potential:'
999 write (iout,'(a/)') 'The epsilon array:'
1000 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1001 write (iout,'(/a)') 'One-body parameters:'
1002 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1003 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1007 C---------------------- GB or BP potential -----------------------------
1009 read (isidep,*,end=116,err=116)(eps(i,j),j=i,ntyp)
1011 read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
1012 read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
1013 read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
1014 read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
1016 c 30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1017 c & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
1018 c & (alp(i),i=1,ntyp)
1019 C For the GB potential convert sigma'**2 into chi'
1022 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1026 write (iout,'(/a/)') 'Parameters of the BP potential:'
1027 write (iout,'(a/)') 'The epsilon array:'
1028 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1029 write (iout,'(/a)') 'One-body parameters:'
1030 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1032 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1033 & chip(i),alp(i),i=1,ntyp)
1036 C--------------------- GBV potential -----------------------------------
1037 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1038 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1039 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1041 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1042 write (iout,'(a/)') 'The epsilon array:'
1043 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1044 write (iout,'(/a)') 'One-body parameters:'
1045 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1046 & 's||/s_|_^2',' chip ',' alph '
1047 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1048 & sigii(i),chip(i),alp(i),i=1,ntyp)
1052 C-----------------------------------------------------------------------
1053 C Calculate the "working" parameters of SC interactions.
1061 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1062 sigma(j,i)=sigma(i,j)
1063 rs0(i,j)=dwa16*sigma(i,j)
1067 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1068 & 'Working parameters of the SC interactions:',
1069 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1074 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1083 sigeps=dsign(1.0D0,epsij)
1085 aa(i,j)=epsij*rrij*rrij
1086 bb(i,j)=-sigeps*epsij*rrij
1090 sigt1sq=sigma0(i)**2
1091 sigt2sq=sigma0(j)**2
1094 ratsig1=sigt2sq/sigt1sq
1095 ratsig2=1.0D0/ratsig1
1096 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1097 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1098 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1102 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1103 sigmaii(i,j)=rsum_max
1104 sigmaii(j,i)=rsum_max
1106 c sigmaii(i,j)=r0(i,j)
1107 c sigmaii(j,i)=r0(i,j)
1109 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1110 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1111 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1112 augm(i,j)=epsij*r_augm**(2*expon)
1113 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1120 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1121 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1122 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1128 C Define the SC-p interaction constants (hard-coded; old style)
1131 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1133 c aad(i,1)=0.3D0*4.0D0**12
1134 C Following line for constants currently implemented
1135 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1136 aad(i,1)=1.5D0*4.0D0**12
1137 c aad(i,1)=0.17D0*5.6D0**12
1139 C "Soft" SC-p repulsion
1141 C Following line for constants currently implemented
1142 c aad(i,1)=0.3D0*4.0D0**6
1143 C "Hard" SC-p repulsion
1144 bad(i,1)=3.0D0*4.0D0**6
1145 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1154 C 8/9/01 Read the SC-p interaction constants from file
1157 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1160 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1161 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1162 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1163 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1167 write (iout,*) "Parameters of SC-p interactions:"
1169 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1170 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1175 C Define the constants of the disulfide bridge
1179 c Old arbitrary potential - commented out.
1184 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1185 c energy surface of diethyl disulfide.
1186 c A. Liwo and U. Kozlowska, 11/24/03
1203 write (iout,'(/a)') "Disulfide bridge parameters:"
1204 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1205 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1206 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1207 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1211 111 write (iout,*) "Error reading bending energy parameters."
1213 112 write (iout,*) "Error reading rotamer energy parameters."
1215 113 write (iout,*) "Error reading torsional energy parameters."
1218 & "Error reading side-chain torsional energy parameters."
1220 114 write (iout,*) "Error reading double torsional energy parameters."
1223 & "Error reading cumulant (multibody energy) parameters."
1225 116 write (iout,*) "Error reading electrostatic energy parameters."
1227 117 write (iout,*) "Error reading side chain interaction parameters."
1229 118 write (iout,*) "Error reading SCp interaction parameters."
1231 119 write (iout,*) "Error reading SCCOR parameters"
1234 call MPI_Finalize(Ierror)
1241 subroutine getenv_loc(var, val)
1242 character(*) var, val
1245 character(2000) line
1248 open (196,file='env',status='old',readonly,shared)
1250 c write(*,*)'looking for ',var
1251 10 read(196,*,err=11,end=11)line
1252 iread=index(line,var)
1253 c write(*,*)iread,' ',var,' ',line
1254 if (iread.eq.0) go to 10
1255 c write(*,*)'---> ',line
1261 iread=iread+ilen(var)+1
1262 read (line(iread:),*,err=12,end=12) val
1263 c write(*,*)'OK: ',var,' = ',val
1269 #elif (defined CRAY)
1270 integer lennam,lenval,ierror
1272 c getenv using a POSIX call, useful on the T3D
1273 c Sept 1996, comment out error check on advice of H. Pritchard
1276 if(lennam.le.0) stop '--error calling getenv--'
1277 call pxfgetenv(var,lennam,val,lenval,ierror)
1278 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1280 call getenv(var,val)