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
248 do j=-nthetyp,nthetyp
249 do k=-nthetyp,nthetyp
250 read (ithep,'(6a)',end=111,err=111) res1
251 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
252 read (ithep,*,end=111,err=111)
253 &(aathet(l,i,j,k,iblock),l=1,ntheterm)
254 read (ithep,*,end=111,err=111)
255 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
256 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
257 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
258 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),
260 read (ithep,*,end=111,err=111)
261 & (((ffthet(llll,lll,ll,i,j,k,iblock),
262 & ffthet(lll,llll,ll,i,j,k,iblock),
263 & ggthet(llll,lll,ll,i,j,k,iblock),
264 & ggthet(lll,llll,ll,i,j,k,iblock),
265 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
270 C For dummy ends assign glycine-type coefficients of theta-only terms; the
271 C coefficients of theta-and-gamma-dependent terms are zero.
276 c aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
277 c aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
279 c aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
280 c aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
283 c aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
285 c aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
291 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
292 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
294 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
295 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
298 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
300 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
304 C Substitution for D aminoacids from symmetry.
307 do j=-nthetyp,nthetyp
308 do k=-nthetyp,nthetyp
309 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
311 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
315 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
316 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
317 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
318 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
324 ffthet(llll,lll,ll,i,j,k,iblock)=
325 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
326 ffthet(lll,llll,ll,i,j,k,iblock)=
327 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
328 ggthet(llll,lll,ll,i,j,k,iblock)=
329 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
330 ggthet(lll,llll,ll,i,j,k,iblock)=
331 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
340 C Control printout of the coefficients of virtual-bond-angle potentials
343 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
347 write (iout,'(//4a)')
348 & 'Type ',onelett(i),onelett(j),onelett(k)
349 write (iout,'(//a,10x,a)') " l","a[l]"
350 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
351 write (iout,'(i2,1pe15.5)')
352 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
354 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))')
355 & "b",l,"c",l,"d",l,"e",l
357 write (iout,'(i2,4(1pe15.5))') m,
358 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
359 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
363 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
364 & "f+",l,"f-",l,"g+",l,"g-",l
367 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
368 & ffthet(n,m,l,i,j,k,iblock),
369 & ffthet(m,n,l,i,j,k,iblock),
370 & ggthet(n,m,l,i,j,k,iblock),
371 & ggthet(m,n,l,i,j,k,iblock)
384 C Read the parameters of the probability distribution/energy expression
385 C of the side chains.
388 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
392 dsc_inv(i)=1.0D0/dsc(i)
403 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
404 & ((blower(k,l,1),l=1,k),k=1,3)
405 censc(1,1,-i)=censc(1,1,i)
406 censc(2,1,-i)=censc(2,1,i)
407 censc(3,1,-i)=-censc(3,1,i)
409 read (irotam,*,end=112,err=112) bsc(j,i)
410 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
411 & ((blower(k,l,j),l=1,k),k=1,3)
412 censc(1,j,-i)=censc(1,j,i)
413 censc(2,j,-i)=censc(2,j,i)
414 censc(3,j,-i)=-censc(3,j,i)
415 C BSC is amplitude of Gaussian
422 akl=akl+blower(k,m,j)*blower(l,m,j)
426 if (((k.eq.3).and.(l.ne.3))
427 & .or.((l.eq.3).and.(k.ne.3))) then
428 gaussc(k,l,j,-i)=-akl
429 gaussc(l,k,j,-i)=-akl
441 write (iout,'(/a)') 'Parameters of side-chain local geometry'
446 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
447 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
448 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
449 & 'log h',(bsc(j,i),j=1,nlobi)
450 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
451 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
453 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
454 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
457 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
458 write (iout,'(a,f10.4,4(16x,f10.4))')
459 & 'Center ',(bsc(j,i),j=1,nlobi)
460 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
469 C Read scrot parameters for potentials determined from all-atom AM1 calculations
470 C added by Urszula Kozlowska 07/11/2007
473 CC TU JEST ZLE musibyc ntyp
474 read (irotam,*,end=112,err=112)
476 read (irotam,*,end=112,err=112)
479 read(irotam,*,end=112,err=112) sc_parmin(j,i)
488 C Read torsional parameters in old format
490 read (itorp,*,end=113,err=113) ntortyp,nterm_old
491 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
492 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
497 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
503 write (iout,'(/a/)') 'Torsional constants:'
506 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
507 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
513 C Read torsional parameters
515 read (itorp,*,end=113,err=113) ntortyp
516 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
519 itortyp(i)=-itortyp(-i)
521 write (iout,*) 'ntortyp',ntortyp
523 do j=-ntortyp+1,ntortyp-1
524 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
526 nterm(-i,-j,iblock)=nterm(i,j,iblock)
527 nlor(-i,-j,iblock)=nlor(i,j,iblock)
530 do k=1,nterm(i,j,iblock)
531 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
533 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
534 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
535 v0ij=v0ij+si*v1(k,i,j,iblock)
537 c write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
538 c write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
539 c &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
541 do k=1,nlor(i,j,iblock)
542 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
543 & vlor2(k,i,j),vlor3(k,i,j)
544 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
547 v0(-i,-j,iblock)=v0ij
553 write (iout,'(/a/)') 'Torsional constants:'
556 write (iout,*) 'ityp',i,' jtyp',j
557 write (iout,*) 'Fourier constants'
558 do k=1,nterm(i,j,iblock)
559 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
562 write (iout,*) 'Lorenz constants'
563 do k=1,nlor(i,j,iblock)
564 write (iout,'(3(1pe15.5))')
565 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
572 C 6/23/01 Read parameters for double torsionals
576 do j=-ntortyp+1,ntortyp-1
577 do k=-ntortyp+1,ntortyp-1
578 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
579 c write (iout,*) "OK onelett",
582 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
583 & .or. t3.ne.toronelet(k)) then
584 write (iout,*) "Error in double torsional parameter file",
587 call MPI_Finalize(Ierror)
589 stop "Error in double torsional parameter file"
591 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
592 & ntermd_2(i,j,k,iblock)
593 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
594 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
595 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
596 & ntermd_1(i,j,k,iblock))
597 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
598 & ntermd_1(i,j,k,iblock))
599 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
600 & ntermd_1(i,j,k,iblock))
601 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
602 & ntermd_1(i,j,k,iblock))
603 C Martix of D parameters for one dimesional foureir series
604 do l=1,ntermd_1(i,j,k,iblock)
605 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
606 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
607 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
608 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
609 c write(iout,*) "whcodze" ,
610 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
612 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
613 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
614 & v2s(m,l,i,j,k,iblock),
615 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
616 C Martix of D parameters for two dimesional fourier series
617 do l=1,ntermd_2(i,j,k,iblock)
619 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
620 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
621 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
622 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
631 write (iout,*) 'Constants for double torsionals'
634 do j=-ntortyp+1,ntortyp-1
635 do k=-ntortyp+1,ntortyp-1
636 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
637 & ' nsingle',ntermd_1(i,j,k,iblock),
638 & ' ndouble',ntermd_2(i,j,k,iblock)
640 write (iout,*) 'Single angles:'
641 do l=1,ntermd_1(i,j,k,iblock)
642 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
643 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
644 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
645 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
648 write (iout,*) 'Pairs of angles:'
649 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
650 do l=1,ntermd_2(i,j,k,iblock)
651 write (iout,'(i5,20f10.5)')
652 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
655 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
656 do l=1,ntermd_2(i,j,k,iblock)
657 write (iout,'(i5,20f10.5)')
658 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
659 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
668 C Read of Side-chain backbone correlation parameters
669 C Modified 11 May 2012 by Adasko
672 read (isccor,*,end=119,err=119) nsccortyp
674 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
676 isccortyp(i)=-isccortyp(-i)
678 iscprol=isccortyp(20)
679 c write (iout,*) 'ntortyp',ntortyp
681 cc maxinter is maximum interaction sites
685 read (isccor,*,end=119,err=119)
686 &nterm_sccor(i,j),nlor_sccor(i,j)
692 nterm_sccor(-i,j)=nterm_sccor(i,j)
693 nterm_sccor(-i,-j)=nterm_sccor(i,j)
694 nterm_sccor(i,-j)=nterm_sccor(i,j)
695 do k=1,nterm_sccor(i,j)
696 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
698 if (j.eq.iscprol) then
699 if (i.eq.isccortyp(10)) then
700 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
701 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
703 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
704 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
705 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
706 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
707 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
708 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
709 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
710 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
713 if (i.eq.isccortyp(10)) then
714 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
715 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
717 if (j.eq.isccortyp(10)) then
718 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
719 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
721 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
722 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
723 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
724 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
725 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
726 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
730 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
731 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
732 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
733 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
736 do k=1,nlor_sccor(i,j)
737 read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
738 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
739 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
740 &(1+vlor3sccor(k,i,j)**2)
742 v0sccor(l,i,j)=v0ijsccor
743 v0sccor(l,-i,j)=v0ijsccor1
744 v0sccor(l,i,-j)=v0ijsccor2
745 v0sccor(l,-i,-j)=v0ijsccor3
751 read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
752 c write (iout,*) 'ntortyp',ntortyp
754 cc maxinter is maximum interaction sites
758 read (isccor,*,end=113,err=113)
759 & nterm_sccor(i,j),nlor_sccor(i,j)
763 do k=1,nterm_sccor(i,j)
764 read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j)
766 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
769 do k=1,nlor_sccor(i,j)
770 read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
771 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
772 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
773 &(1+vlor3sccor(k,i,j)**2)
775 v0sccor(i,j,iblock)=v0ijsccor
783 write (iout,'(/a/)') 'Torsional constants:'
786 write (iout,*) 'ityp',i,' jtyp',j
787 write (iout,*) 'Fourier constants'
788 do k=1,nterm_sccor(i,j)
789 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
791 write (iout,*) 'Lorenz constants'
792 do k=1,nlor_sccor(i,j)
793 write (iout,'(3(1pe15.5))')
794 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
801 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
802 C interaction energy of the Gly, Ala, and Pro prototypes.
806 write (iout,*) "Coefficients of the cumulants"
808 read (ifourier,*) nloctyp
810 read (ifourier,*,end=115,err=115)
811 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
813 write (iout,*) 'Type',i
814 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
856 c Ctilde(1,1,i)=0.0d0
857 c Ctilde(1,2,i)=0.0d0
858 c Ctilde(2,1,i)=0.0d0
859 c Ctilde(2,2,i)=0.0d0
881 c Dtilde(1,1,i)=0.0d0
882 c Dtilde(1,2,i)=0.0d0
883 c Dtilde(2,1,i)=0.0d0
884 c Dtilde(2,2,i)=0.0d0
885 EE(1,1,i)= b(10)+b(11)
886 EE(2,2,i)=-b(10)+b(11)
887 EE(2,1,i)= b(12)-b(13)
888 EE(1,2,i)= b(12)+b(13)
889 EE(1,1,-i)= b(10)+b(11)
890 EE(2,2,-i)=-b(10)+b(11)
891 EE(2,1,-i)=-b(12)+b(13)
892 EE(1,2,-i)=-b(12)-b(13)
898 c ee(2,1,i)=ee(1,2,i)
902 write (iout,*) 'Type',i
904 write(iout,*) B1(1,i),B1(2,i)
906 write(iout,*) B2(1,i),B2(2,i)
909 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
913 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
917 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
923 C Read electrostatic-interaction parameters
927 write (iout,'(/a)') 'Electrostatic interaction constants:'
928 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
929 & 'IT','JT','APP','BPP','AEL6','AEL3'
931 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
932 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
933 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
934 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
939 app (i,j)=epp(i,j)*rri*rri
940 bpp (i,j)=-2.0D0*epp(i,j)*rri
941 ael6(i,j)=elpp6(i,j)*4.2D0**6
942 ael3(i,j)=elpp3(i,j)*4.2D0**3
944 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
945 & ael6(i,j),ael3(i,j)
950 C Read side-chain interaction parameters.
952 read (isidep,*,end=117,err=117) ipot,expon
953 if (ipot.lt.1 .or. ipot.gt.5) then
954 write (iout,'(2a)') 'Error while reading SC interaction',
955 & 'potential file - unknown potential type.'
957 call MPI_Finalize(Ierror)
963 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
964 & ', exponents are ',expon,2*expon
965 goto (10,20,30,30,40) ipot
966 C----------------------- LJ potential ---------------------------------
967 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
968 & (sigma0(i),i=1,ntyp)
970 write (iout,'(/a/)') 'Parameters of the LJ potential:'
971 write (iout,'(a/)') 'The epsilon array:'
972 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
973 write (iout,'(/a)') 'One-body parameters:'
974 write (iout,'(a,4x,a)') 'residue','sigma'
975 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
978 C----------------------- LJK potential --------------------------------
979 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
980 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
982 write (iout,'(/a/)') 'Parameters of the LJK potential:'
983 write (iout,'(a/)') 'The epsilon array:'
984 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
985 write (iout,'(/a)') 'One-body parameters:'
986 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
987 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
991 C---------------------- GB or BP potential -----------------------------
993 read (isidep,*,end=116,err=116)(eps(i,j),j=i,ntyp)
995 read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
996 read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
997 read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
998 read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
999 C For the GB potential convert sigma'**2 into chi'
1002 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1006 write (iout,'(/a/)') 'Parameters of the BP potential:'
1007 write (iout,'(a/)') 'The epsilon array:'
1008 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1009 write (iout,'(/a)') 'One-body parameters:'
1010 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1012 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1013 & chip(i),alp(i),i=1,ntyp)
1016 C--------------------- GBV potential -----------------------------------
1017 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1018 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1019 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1021 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1022 write (iout,'(a/)') 'The epsilon array:'
1023 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1024 write (iout,'(/a)') 'One-body parameters:'
1025 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1026 & 's||/s_|_^2',' chip ',' alph '
1027 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1028 & sigii(i),chip(i),alp(i),i=1,ntyp)
1032 C-----------------------------------------------------------------------
1033 C Calculate the "working" parameters of SC interactions.
1041 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1042 sigma(j,i)=sigma(i,j)
1043 rs0(i,j)=dwa16*sigma(i,j)
1047 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1048 & 'Working parameters of the SC interactions:',
1049 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1054 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1063 sigeps=dsign(1.0D0,epsij)
1065 aa(i,j)=epsij*rrij*rrij
1066 bb(i,j)=-sigeps*epsij*rrij
1070 sigt1sq=sigma0(i)**2
1071 sigt2sq=sigma0(j)**2
1074 ratsig1=sigt2sq/sigt1sq
1075 ratsig2=1.0D0/ratsig1
1076 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1077 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1078 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1082 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1083 sigmaii(i,j)=rsum_max
1084 sigmaii(j,i)=rsum_max
1086 c sigmaii(i,j)=r0(i,j)
1087 c sigmaii(j,i)=r0(i,j)
1089 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1090 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1091 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1092 augm(i,j)=epsij*r_augm**(2*expon)
1093 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1100 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1101 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1102 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1108 C Define the SC-p interaction constants (hard-coded; old style)
1111 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1113 c aad(i,1)=0.3D0*4.0D0**12
1114 C Following line for constants currently implemented
1115 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1116 aad(i,1)=1.5D0*4.0D0**12
1117 c aad(i,1)=0.17D0*5.6D0**12
1119 C "Soft" SC-p repulsion
1121 C Following line for constants currently implemented
1122 c aad(i,1)=0.3D0*4.0D0**6
1123 C "Hard" SC-p repulsion
1124 bad(i,1)=3.0D0*4.0D0**6
1125 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1134 C 8/9/01 Read the SC-p interaction constants from file
1137 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1140 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1141 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1142 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1143 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1147 write (iout,*) "Parameters of SC-p interactions:"
1149 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1150 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1156 C Define the constants of the disulfide bridge
1160 c Old arbitrary potential - commented out.
1165 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1166 c energy surface of diethyl disulfide.
1167 c A. Liwo and U. Kozlowska, 11/24/03
1184 write (iout,'(/a)') "Disulfide bridge parameters:"
1185 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1186 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1187 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1188 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1192 111 write (iout,*) "Error reading bending energy parameters."
1194 112 write (iout,*) "Error reading rotamer energy parameters."
1196 113 write (iout,*) "Error reading torsional energy parameters."
1198 114 write (iout,*) "Error reading double torsional energy parameters."
1201 & "Error reading cumulant (multibody energy) parameters."
1203 116 write (iout,*) "Error reading electrostatic energy parameters."
1205 117 write (iout,*) "Error reading side chain interaction parameters."
1207 118 write (iout,*) "Error reading SCp interaction parameters."
1209 119 write (iout,*) "Error reading SCCOR parameters"
1212 call MPI_Finalize(Ierror)
1219 subroutine getenv_loc(var, val)
1220 character(*) var, val
1223 character(2000) line
1226 open (196,file='env',status='old',readonly,shared)
1228 c write(*,*)'looking for ',var
1229 10 read(196,*,err=11,end=11)line
1230 iread=index(line,var)
1231 c write(*,*)iread,' ',var,' ',line
1232 if (iread.eq.0) go to 10
1233 c write(*,*)'---> ',line
1239 iread=iread+ilen(var)+1
1240 read (line(iread:),*,err=12,end=12) val
1241 c write(*,*)'OK: ',var,' = ',val
1247 #elif (defined CRAY)
1248 integer lennam,lenval,ierror
1250 c getenv using a POSIX call, useful on the T3D
1251 c Sept 1996, comment out error check on advice of H. Pritchard
1254 if(lennam.le.0) stop '--error calling getenv--'
1255 call pxfgetenv(var,lennam,val,lenval,ierror)
1256 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1258 call getenv(var,val)