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))
712 C Read of Side-chain backbone correlation parameters
713 C Modified 11 May 2012 by Adasko
716 read (isccor,*,end=1113,err=1113) nsccortyp
718 read (isccor,*,end=1113,err=1113) (isccortyp(i),i=1,ntyp)
720 isccortyp(i)=-isccortyp(-i)
722 iscprol=isccortyp(20)
723 c write (iout,*) 'ntortyp',ntortyp
725 cc maxinter is maximum interaction sites
729 read (isccor,*,end=1113,err=1113) nterm_sccor(i,j),
736 nterm_sccor(-i,j)=nterm_sccor(i,j)
737 nterm_sccor(-i,-j)=nterm_sccor(i,j)
738 nterm_sccor(i,-j)=nterm_sccor(i,j)
739 do k=1,nterm_sccor(i,j)
740 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
742 if (j.eq.iscprol) then
743 if (i.eq.isccortyp(10)) then
744 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
745 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
747 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
748 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
749 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
750 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
751 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
752 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
753 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
754 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
757 if (i.eq.isccortyp(10)) then
758 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
759 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
761 if (j.eq.isccortyp(10)) then
762 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
763 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
765 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
766 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
767 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
768 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
769 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
770 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
774 read (isccor,*,end=1113,err=1113) kk,v1sccor(k,l,i,j)
776 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
777 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
778 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
779 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
782 do k=1,nlor_sccor(i,j)
783 read (isccor,*,end=1113,err=1113) kk,vlor1sccor(k,i,j),
784 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
785 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
786 &(1+vlor3sccor(k,i,j)**2)
788 v0sccor(l,i,j)=v0ijsccor
789 v0sccor(l,-i,j)=v0ijsccor1
790 v0sccor(l,i,-j)=v0ijsccor2
791 v0sccor(l,-i,-j)=v0ijsccor3
797 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
798 c write (iout,*) 'ntortyp',ntortyp
800 cc maxinter is maximum interaction sites
804 read (isccor,*,end=119,err=119)
805 & nterm_sccor(i,j),nlor_sccor(i,j)
809 do k=1,nterm_sccor(i,j)
810 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
812 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
815 do k=1,nlor_sccor(i,j)
816 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
817 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
818 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
819 &(1+vlor3sccor(k,i,j)**2)
821 v0sccor(i,j,iblock)=v0ijsccor
829 write (iout,'(/a/)') 'Torsional constants:'
832 write (iout,*) 'ityp',i,' jtyp',j
833 write (iout,*) 'Fourier constants'
834 do k=1,nterm_sccor(i,j)
835 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
837 write (iout,*) 'Lorenz constants'
838 do k=1,nlor_sccor(i,j)
839 write (iout,'(3(1pe15.5))')
840 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
847 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
848 C interaction energy of the Gly, Ala, and Pro prototypes.
852 write (iout,*) "Coefficients of the cumulants"
854 read (ifourier,*) nloctyp
856 read (ifourier,*,end=115,err=115)
857 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
859 write (iout,*) 'Type',i
860 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
888 c Ctilde(1,1,i)=0.0d0
889 c Ctilde(1,2,i)=0.0d0
890 c Ctilde(2,1,i)=0.0d0
891 c Ctilde(2,2,i)=0.0d0
904 c Dtilde(1,1,i)=0.0d0
905 c Dtilde(1,2,i)=0.0d0
906 c Dtilde(2,1,i)=0.0d0
907 c Dtilde(2,2,i)=0.0d0
908 EE(1,1,i)= b(10)+b(11)
909 EE(2,2,i)=-b(10)+b(11)
910 EE(2,1,i)= b(12)-b(13)
911 EE(1,2,i)= b(12)+b(13)
916 c ee(2,1,i)=ee(1,2,i)
920 write (iout,*) 'Type',i
922 write(iout,*) B1(1,i),B1(2,i)
924 write(iout,*) B2(1,i),B2(2,i)
927 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
931 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
935 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
940 C Read electrostatic-interaction parameters
944 write (iout,'(/a)') 'Electrostatic interaction constants:'
945 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
946 & 'IT','JT','APP','BPP','AEL6','AEL3'
948 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
949 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
950 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
951 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
956 app (i,j)=epp(i,j)*rri*rri
957 bpp (i,j)=-2.0D0*epp(i,j)*rri
958 ael6(i,j)=elpp6(i,j)*4.2D0**6
959 ael3(i,j)=elpp3(i,j)*4.2D0**3
960 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
961 & ael6(i,j),ael3(i,j)
965 C Read side-chain interaction parameters.
967 read (isidep,*,end=117,err=117) ipot,expon
968 if (ipot.lt.1 .or. ipot.gt.5) then
969 write (iout,'(2a)') 'Error while reading SC interaction',
970 & 'potential file - unknown potential type.'
972 call MPI_Finalize(Ierror)
978 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
979 & ', exponents are ',expon,2*expon
980 goto (10,20,30,30,40) ipot
981 C----------------------- LJ potential ---------------------------------
982 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
983 & (sigma0(i),i=1,ntyp)
985 write (iout,'(/a/)') 'Parameters of the LJ potential:'
986 write (iout,'(a/)') 'The epsilon array:'
987 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
988 write (iout,'(/a)') 'One-body parameters:'
989 write (iout,'(a,4x,a)') 'residue','sigma'
990 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
993 C----------------------- LJK potential --------------------------------
994 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
995 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
997 write (iout,'(/a/)') 'Parameters of the LJK potential:'
998 write (iout,'(a/)') 'The epsilon array:'
999 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1000 write (iout,'(/a)') 'One-body parameters:'
1001 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1002 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1006 C---------------------- GB or BP potential -----------------------------
1008 read (isidep,*,end=116,err=116)(eps(i,j),j=i,ntyp)
1010 read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
1011 read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
1012 read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
1013 read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
1015 c 30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1016 c & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
1017 c & (alp(i),i=1,ntyp)
1018 C For the GB potential convert sigma'**2 into chi'
1021 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1025 write (iout,'(/a/)') 'Parameters of the BP potential:'
1026 write (iout,'(a/)') 'The epsilon array:'
1027 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1028 write (iout,'(/a)') 'One-body parameters:'
1029 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1031 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1032 & chip(i),alp(i),i=1,ntyp)
1035 C--------------------- GBV potential -----------------------------------
1036 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1037 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1038 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1040 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1041 write (iout,'(a/)') 'The epsilon array:'
1042 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1043 write (iout,'(/a)') 'One-body parameters:'
1044 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1045 & 's||/s_|_^2',' chip ',' alph '
1046 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1047 & sigii(i),chip(i),alp(i),i=1,ntyp)
1051 C-----------------------------------------------------------------------
1052 C Calculate the "working" parameters of SC interactions.
1060 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1061 sigma(j,i)=sigma(i,j)
1062 rs0(i,j)=dwa16*sigma(i,j)
1066 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1067 & 'Working parameters of the SC interactions:',
1068 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1073 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1082 sigeps=dsign(1.0D0,epsij)
1084 aa(i,j)=epsij*rrij*rrij
1085 bb(i,j)=-sigeps*epsij*rrij
1089 sigt1sq=sigma0(i)**2
1090 sigt2sq=sigma0(j)**2
1093 ratsig1=sigt2sq/sigt1sq
1094 ratsig2=1.0D0/ratsig1
1095 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1096 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1097 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1101 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1102 sigmaii(i,j)=rsum_max
1103 sigmaii(j,i)=rsum_max
1105 c sigmaii(i,j)=r0(i,j)
1106 c sigmaii(j,i)=r0(i,j)
1108 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1109 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1110 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1111 augm(i,j)=epsij*r_augm**(2*expon)
1112 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1119 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1120 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1121 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1127 C Define the SC-p interaction constants (hard-coded; old style)
1130 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1132 c aad(i,1)=0.3D0*4.0D0**12
1133 C Following line for constants currently implemented
1134 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1135 aad(i,1)=1.5D0*4.0D0**12
1136 c aad(i,1)=0.17D0*5.6D0**12
1138 C "Soft" SC-p repulsion
1140 C Following line for constants currently implemented
1141 c aad(i,1)=0.3D0*4.0D0**6
1142 C "Hard" SC-p repulsion
1143 bad(i,1)=3.0D0*4.0D0**6
1144 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1153 C 8/9/01 Read the SC-p interaction constants from file
1156 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1159 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1160 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1161 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1162 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1166 write (iout,*) "Parameters of SC-p interactions:"
1168 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1169 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1174 C Define the constants of the disulfide bridge
1178 c Old arbitrary potential - commented out.
1183 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1184 c energy surface of diethyl disulfide.
1185 c A. Liwo and U. Kozlowska, 11/24/03
1202 write (iout,'(/a)') "Disulfide bridge parameters:"
1203 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1204 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1205 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1206 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1210 111 write (iout,*) "Error reading bending energy parameters."
1212 112 write (iout,*) "Error reading rotamer energy parameters."
1214 113 write (iout,*) "Error reading torsional energy parameters."
1217 & "Error reading side-chain torsional energy parameters."
1219 114 write (iout,*) "Error reading double torsional energy parameters."
1222 & "Error reading cumulant (multibody energy) parameters."
1224 116 write (iout,*) "Error reading electrostatic energy parameters."
1226 117 write (iout,*) "Error reading side chain interaction parameters."
1228 118 write (iout,*) "Error reading SCp interaction parameters."
1230 119 write (iout,*) "Error reading SCCOR parameters"
1233 call MPI_Finalize(Ierror)
1240 subroutine getenv_loc(var, val)
1241 character(*) var, val
1244 character(2000) line
1247 open (196,file='env',status='old',readonly,shared)
1249 c write(*,*)'looking for ',var
1250 10 read(196,*,err=11,end=11)line
1251 iread=index(line,var)
1252 c write(*,*)iread,' ',var,' ',line
1253 if (iread.eq.0) go to 10
1254 c write(*,*)'---> ',line
1260 iread=iread+ilen(var)+1
1261 read (line(iread:),*,err=12,end=12) val
1262 c write(*,*)'OK: ',var,' = ',val
1268 #elif (defined CRAY)
1269 integer lennam,lenval,ierror
1271 c getenv using a POSIX call, useful on the T3D
1272 c Sept 1996, comment out error check on advice of H. Pritchard
1275 if(lennam.le.0) stop '--error calling getenv--'
1276 call pxfgetenv(var,lennam,val,lenval,ierror)
1277 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1279 call getenv(var,val)