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
214 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
215 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
217 ithetyp(i)=-ithetyp(-i)
220 do i=-maxthetyp,maxthetyp
221 do j=-maxthetyp,maxthetyp
222 do k=-maxthetyp,maxthetyp
223 aa0thet(i,j,k,iblock)=0.0d0
225 aathet(l,i,j,k,iblock)=0.0d0
229 bbthet(m,l,i,j,k,iblock)=0.0d0
230 ccthet(m,l,i,j,k,iblock)=0.0d0
231 ddthet(m,l,i,j,k,iblock)=0.0d0
232 eethet(m,l,i,j,k,iblock)=0.0d0
238 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
239 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
247 c VAR:iblock means terminally blocking group 1=non-proline 2=proline
249 c VAR:ntethtyp is type of theta potentials type currently 0=glycine
250 c VAR:1=non-glicyne non-proline 2=proline
251 c VAR:negative values for D-aminoacid
254 do j=-nthetyp,nthetyp
255 do k=-nthetyp,nthetyp
256 read (ithep,'(6a)',end=111,err=111) res1
257 c VAR: aa0thet is variable describing the average value of Foureir
258 c VAR: expansion series
259 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
260 c VAR: aathet is foureir expansion in theta/2 angle for full formula
261 c VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
262 Condens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
263 read (ithep,*,end=111,err=111)
264 &(aathet(l,i,j,k,iblock),l=1,ntheterm)
265 read (ithep,*,end=111,err=111)
266 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
267 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
268 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
269 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),
271 read (ithep,*,end=111,err=111)
272 & (((ffthet(llll,lll,ll,i,j,k,iblock),
273 & ffthet(lll,llll,ll,i,j,k,iblock),
274 & ggthet(llll,lll,ll,i,j,k,iblock),
275 & ggthet(lll,llll,ll,i,j,k,iblock),
276 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
283 C For dummy ends assign glycine-type coefficients of theta-only terms; the
284 C coefficients of theta-and-gamma-dependent terms are zero.
285 C IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
286 C RECOMENTDED AFTER VERSION 3.3)
290 c aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
291 c aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
293 c aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
294 c aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
297 c aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
299 c aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
302 C AND COMMENT THE LOOPS BELOW
306 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
307 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
309 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
310 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
313 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
315 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
319 C Substitution for D aminoacids from symmetry.
322 do j=-nthetyp,nthetyp
323 do k=-nthetyp,nthetyp
324 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
326 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
330 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
331 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
332 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
333 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
339 ffthet(llll,lll,ll,i,j,k,iblock)=
340 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
341 ffthet(lll,llll,ll,i,j,k,iblock)=
342 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
343 ggthet(llll,lll,ll,i,j,k,iblock)=
344 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
345 ggthet(lll,llll,ll,i,j,k,iblock)=
346 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
355 C Control printout of the coefficients of virtual-bond-angle potentials
358 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
362 write (iout,'(//4a)')
363 & 'Type ',onelett(i),onelett(j),onelett(k)
364 write (iout,'(//a,10x,a)') " l","a[l]"
365 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
366 write (iout,'(i2,1pe15.5)')
367 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
369 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))')
370 & "b",l,"c",l,"d",l,"e",l
372 write (iout,'(i2,4(1pe15.5))') m,
373 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
374 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
378 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
379 & "f+",l,"f-",l,"g+",l,"g-",l
382 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
383 & ffthet(n,m,l,i,j,k,iblock),
384 & ffthet(m,n,l,i,j,k,iblock),
385 & ggthet(n,m,l,i,j,k,iblock),
386 & ggthet(m,n,l,i,j,k,iblock)
395 write (2,*) "Start reading THETA_PDB"
397 read (ithep_pdb,*,err=111,end=111)
398 & a0thet(i),(athet(j,i,1,1),j=1,2),
399 & (bthet(j,i,1,1),j=1,2)
400 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
401 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
402 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
406 athet(1,i,1,-1)=athet(1,i,1,1)
407 athet(2,i,1,-1)=athet(2,i,1,1)
408 bthet(1,i,1,-1)=-bthet(1,i,1,1)
409 bthet(2,i,1,-1)=-bthet(2,i,1,1)
410 athet(1,i,-1,1)=-athet(1,i,1,1)
411 athet(2,i,-1,1)=-athet(2,i,1,1)
412 bthet(1,i,-1,1)=bthet(1,i,1,1)
413 bthet(2,i,-1,1)=bthet(2,i,1,1)
417 athet(1,i,-1,-1)=athet(1,-i,1,1)
418 athet(2,i,-1,-1)=-athet(2,-i,1,1)
419 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
420 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
421 athet(1,i,-1,1)=athet(1,-i,1,1)
422 athet(2,i,-1,1)=-athet(2,-i,1,1)
423 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
424 bthet(2,i,-1,1)=bthet(2,-i,1,1)
425 athet(1,i,1,-1)=-athet(1,-i,1,1)
426 athet(2,i,1,-1)=athet(2,-i,1,1)
427 bthet(1,i,1,-1)=bthet(1,-i,1,1)
428 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
433 polthet(j,i)=polthet(j,-i)
436 gthet(j,i)=gthet(j,-i)
439 write (2,*) "End reading THETA_PDB"
445 C Read the parameters of the probability distribution/energy expression
446 C of the side chains.
449 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
453 dsc_inv(i)=1.0D0/dsc(i)
464 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
465 & ((blower(k,l,1),l=1,k),k=1,3)
467 read (irotam,*,end=112,err=112) bsc(j,i)
468 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
469 & ((blower(k,l,j),l=1,k),k=1,3)
476 akl=akl+blower(k,m,j)*blower(l,m,j)
487 write (iout,'(/a)') 'Parameters of side-chain local geometry'
492 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
493 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
494 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
495 & 'log h',(bsc(j,i),j=1,nlobi)
496 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
497 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
499 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
500 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
503 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
504 write (iout,'(a,f10.4,4(16x,f10.4))')
505 & 'Center ',(bsc(j,i),j=1,nlobi)
506 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
515 C Read scrot parameters for potentials determined from all-atom AM1 calculations
516 C added by Urszula Kozlowska 07/11/2007
519 read (irotam,*,end=112,err=112)
521 read (irotam,*,end=112,err=112)
524 read(irotam,*,end=112,err=112) sc_parmin(j,i)
529 C Read the parameters of the probability distribution/energy expression
530 C of the side chains.
532 write (2,*) "Start reading ROTAM_PDB"
535 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
539 dsc_inv(i)=1.0D0/dsc(i)
550 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
551 & ((blower(k,l,1),l=1,k),k=1,3)
553 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
554 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
555 & ((blower(k,l,j),l=1,k),k=1,3)
562 akl=akl+blower(k,m,j)*blower(l,m,j)
572 write (2,*) "Ending reading ROTAM_PDB"
578 C Read torsional parameters in old format
580 read (itorp,*,end=113,err=113) ntortyp,nterm_old
581 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
582 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
587 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
593 write (iout,'(/a/)') 'Torsional constants:'
596 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
597 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
603 C Read torsional parameters
605 read (itorp,*,end=113,err=113) ntortyp
606 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
607 c write (iout,*) 'ntortyp',ntortyp
610 read (itorp,*,end=113,err=113) nterm(i,j),nlor(i,j)
614 read (itorp,*,end=113,err=113) kk,v1(k,i,j),v2(k,i,j)
615 v0ij=v0ij+si*v1(k,i,j)
619 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
620 & vlor2(k,i,j),vlor3(k,i,j)
621 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
628 write (iout,'(/a/)') 'Torsional constants:'
631 write (iout,*) 'ityp',i,' jtyp',j
632 write (iout,*) 'Fourier constants'
634 write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
636 write (iout,*) 'Lorenz constants'
638 write (iout,'(3(1pe15.5))')
639 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
645 C 6/23/01 Read parameters for double torsionals
650 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
651 write (iout,*) "OK onelett",
654 if (t1.ne.onelett(i) .or. t2.ne.onelett(j)
655 & .or. t3.ne.onelett(k)) then
656 write (iout,*) "Error in double torsional parameter file",
659 call MPI_Finalize(Ierror)
661 stop "Error in double torsional parameter file"
663 read (itordp,*,end=114,err=114) ntermd_1(i,j,k),
665 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k),l=1,
667 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k),l=1,
669 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k),l=1,
671 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k),l=1,
673 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k),
674 & v2c(m,l,i,j,k),v2s(l,m,i,j,k),v2s(m,l,i,j,k),
675 & m=1,l-1),l=1,ntermd_2(i,j,k))
681 write (iout,*) 'Constants for double torsionals'
684 do j=-ntortyp+1,ntortyp-1
685 do k=-ntortyp+1,ntortyp-1
686 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
687 & ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k)
689 write (iout,*) 'Single angles:'
690 do l=1,ntermd_1(i,j,k)
691 write (iout,'(i5,2f10.5,5x,2f10.5)') l,
692 & v1c(1,l,i,j,k),v1s(1,l,i,j,k),
693 & v1c(2,l,i,j,k),v1s(2,l,i,j,k)
696 write (iout,*) 'Pairs of angles:'
697 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
698 do l=1,ntermd_2(i,j,k)
699 write (iout,'(i5,20f10.5)')
700 & l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k))
703 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
704 do l=1,ntermd_2(i,j,k)
705 write (iout,'(i5,20f10.5)')
706 & l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k))
715 C Read of Side-chain backbone correlation parameters
716 C Modified 11 May 2012 by Adasko
719 read (isccor,*,end=1113,err=1113) nsccortyp
721 read (isccor,*,end=1113,err=1113) (isccortyp(i),i=1,ntyp)
723 isccortyp(i)=-isccortyp(-i)
725 iscprol=isccortyp(20)
726 c write (iout,*) 'ntortyp',ntortyp
728 cc maxinter is maximum interaction sites
732 read (isccor,*,end=1113,err=1113) nterm_sccor(i,j),
740 nterm_sccor(-i,j)=nterm_sccor(i,j)
741 nterm_sccor(-i,-j)=nterm_sccor(i,j)
742 nterm_sccor(i,-j)=nterm_sccor(i,j)
743 do k=1,nterm_sccor(i,j)
744 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
746 if (j.eq.iscprol) then
747 if (i.eq.isccortyp(10)) then
748 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
749 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
751 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
752 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
753 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
754 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
755 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
756 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
757 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
758 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
761 if (i.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 if (j.eq.isccortyp(10)) then
766 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
767 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)
771 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
772 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
773 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
774 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
778 C read (isccor,*,end=1113,err=1113) kk,v1sccor(k,l,i,j)
779 C & ,v2sccor(k,l,i,j)
780 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
781 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
782 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
783 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
786 do k=1,nlor_sccor(i,j)
787 read (isccor,*,end=1113,err=1113) kk,vlor1sccor(k,i,j),
788 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
789 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
790 &(1+vlor3sccor(k,i,j)**2)
792 v0sccor(l,i,j)=v0ijsccor
793 v0sccor(l,-i,j)=v0ijsccor1
794 v0sccor(l,i,-j)=v0ijsccor2
795 v0sccor(l,-i,-j)=v0ijsccor3
801 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
802 c write (iout,*) 'ntortyp',ntortyp
804 cc maxinter is maximum interaction sites
808 read (isccor,*,end=119,err=119)
809 & nterm_sccor(i,j),nlor_sccor(i,j)
813 do k=1,nterm_sccor(i,j)
814 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
816 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
819 do k=1,nlor_sccor(i,j)
820 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
821 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
822 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
823 &(1+vlor3sccor(k,i,j)**2)
825 v0sccor(i,j,iblock)=v0ijsccor
833 write (iout,'(/a/)') 'Torsional constants:'
836 write (iout,*) 'ityp',i,' jtyp',j
837 write (iout,*) 'Fourier constants'
838 do k=1,nterm_sccor(i,j)
839 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
841 write (iout,*) 'Lorenz constants'
842 do k=1,nlor_sccor(i,j)
843 write (iout,'(3(1pe15.5))')
844 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
851 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
852 C interaction energy of the Gly, Ala, and Pro prototypes.
856 write (iout,*) "Coefficients of the cumulants"
858 read (ifourier,*) nloctyp
860 read (ifourier,*,end=115,err=115)
861 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
863 write (iout,*) 'Type',i
864 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
872 C B1tilde(1,-i) =-b(3)
873 C B1tilde(2,-i) =b(5)
892 c Ctilde(1,1,i)=0.0d0
893 c Ctilde(1,2,i)=0.0d0
894 c Ctilde(2,1,i)=0.0d0
895 c Ctilde(2,2,i)=0.0d0
908 c Dtilde(1,1,i)=0.0d0
909 c Dtilde(1,2,i)=0.0d0
910 c Dtilde(2,1,i)=0.0d0
911 c Dtilde(2,2,i)=0.0d0
912 EE(1,1,i)= b(10)+b(11)
913 EE(2,2,i)=-b(10)+b(11)
914 EE(2,1,i)= b(12)-b(13)
915 EE(1,2,i)= b(12)+b(13)
920 c ee(2,1,i)=ee(1,2,i)
924 write (iout,*) 'Type',i
926 write(iout,*) B1(1,i),B1(2,i)
928 write(iout,*) B2(1,i),B2(2,i)
931 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
935 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
939 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
944 C Read electrostatic-interaction parameters
948 write (iout,'(/a)') 'Electrostatic interaction constants:'
949 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
950 & 'IT','JT','APP','BPP','AEL6','AEL3'
952 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
953 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
954 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
955 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
960 app (i,j)=epp(i,j)*rri*rri
961 bpp (i,j)=-2.0D0*epp(i,j)*rri
962 ael6(i,j)=elpp6(i,j)*4.2D0**6
963 ael3(i,j)=elpp3(i,j)*4.2D0**3
964 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
965 & ael6(i,j),ael3(i,j)
969 C Read side-chain interaction parameters.
971 read (isidep,*,end=117,err=117) ipot,expon
972 if (ipot.lt.1 .or. ipot.gt.5) then
973 write (iout,'(2a)') 'Error while reading SC interaction',
974 & 'potential file - unknown potential type.'
976 call MPI_Finalize(Ierror)
982 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
983 & ', exponents are ',expon,2*expon
984 goto (10,20,30,30,40) ipot
985 C----------------------- LJ potential ---------------------------------
986 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
987 & (sigma0(i),i=1,ntyp)
989 write (iout,'(/a/)') 'Parameters of the LJ potential:'
990 write (iout,'(a/)') 'The epsilon array:'
991 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
992 write (iout,'(/a)') 'One-body parameters:'
993 write (iout,'(a,4x,a)') 'residue','sigma'
994 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
997 C----------------------- LJK potential --------------------------------
998 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
999 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1001 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1002 write (iout,'(a/)') 'The epsilon array:'
1003 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1004 write (iout,'(/a)') 'One-body parameters:'
1005 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1006 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1010 C---------------------- GB or BP potential -----------------------------
1012 read (isidep,*,end=116,err=116)(eps(i,j),j=i,ntyp)
1014 read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
1015 read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
1016 read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
1017 read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
1019 c 30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1020 c & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
1021 c & (alp(i),i=1,ntyp)
1022 C For the GB potential convert sigma'**2 into chi'
1025 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1029 write (iout,'(/a/)') 'Parameters of the BP potential:'
1030 write (iout,'(a/)') 'The epsilon array:'
1031 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1032 write (iout,'(/a)') 'One-body parameters:'
1033 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1035 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1036 & chip(i),alp(i),i=1,ntyp)
1039 C--------------------- GBV potential -----------------------------------
1040 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1041 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1042 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1044 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1045 write (iout,'(a/)') 'The epsilon array:'
1046 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1047 write (iout,'(/a)') 'One-body parameters:'
1048 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1049 & 's||/s_|_^2',' chip ',' alph '
1050 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1051 & sigii(i),chip(i),alp(i),i=1,ntyp)
1055 C-----------------------------------------------------------------------
1056 C Calculate the "working" parameters of SC interactions.
1064 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1065 sigma(j,i)=sigma(i,j)
1066 rs0(i,j)=dwa16*sigma(i,j)
1070 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1071 & 'Working parameters of the SC interactions:',
1072 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1077 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1086 sigeps=dsign(1.0D0,epsij)
1088 aa(i,j)=epsij*rrij*rrij
1089 bb(i,j)=-sigeps*epsij*rrij
1093 sigt1sq=sigma0(i)**2
1094 sigt2sq=sigma0(j)**2
1097 ratsig1=sigt2sq/sigt1sq
1098 ratsig2=1.0D0/ratsig1
1099 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1100 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1101 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1105 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1106 sigmaii(i,j)=rsum_max
1107 sigmaii(j,i)=rsum_max
1109 c sigmaii(i,j)=r0(i,j)
1110 c sigmaii(j,i)=r0(i,j)
1112 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1113 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1114 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1115 augm(i,j)=epsij*r_augm**(2*expon)
1116 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1123 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1124 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1125 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1131 C Define the SC-p interaction constants (hard-coded; old style)
1134 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1136 c aad(i,1)=0.3D0*4.0D0**12
1137 C Following line for constants currently implemented
1138 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1139 aad(i,1)=1.5D0*4.0D0**12
1140 c aad(i,1)=0.17D0*5.6D0**12
1142 C "Soft" SC-p repulsion
1144 C Following line for constants currently implemented
1145 c aad(i,1)=0.3D0*4.0D0**6
1146 C "Hard" SC-p repulsion
1147 bad(i,1)=3.0D0*4.0D0**6
1148 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1157 C 8/9/01 Read the SC-p interaction constants from file
1160 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1163 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1164 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1165 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1166 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1170 write (iout,*) "Parameters of SC-p interactions:"
1172 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1173 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1178 C Define the constants of the disulfide bridge
1182 c Old arbitrary potential - commented out.
1187 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1188 c energy surface of diethyl disulfide.
1189 c A. Liwo and U. Kozlowska, 11/24/03
1206 write (iout,'(/a)') "Disulfide bridge parameters:"
1207 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1208 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1209 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1210 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1214 111 write (iout,*) "Error reading bending energy parameters."
1216 112 write (iout,*) "Error reading rotamer energy parameters."
1218 113 write (iout,*) "Error reading torsional energy parameters."
1221 & "Error reading side-chain torsional energy parameters."
1223 114 write (iout,*) "Error reading double torsional energy parameters."
1226 & "Error reading cumulant (multibody energy) parameters."
1228 116 write (iout,*) "Error reading electrostatic energy parameters."
1230 117 write (iout,*) "Error reading side chain interaction parameters."
1232 118 write (iout,*) "Error reading SCp interaction parameters."
1234 119 write (iout,*) "Error reading SCCOR parameters"
1237 call MPI_Finalize(Ierror)
1244 subroutine getenv_loc(var, val)
1245 character(*) var, val
1248 character(2000) line
1251 open (196,file='env',status='old',readonly,shared)
1253 c write(*,*)'looking for ',var
1254 10 read(196,*,err=11,end=11)line
1255 iread=index(line,var)
1256 c write(*,*)iread,' ',var,' ',line
1257 if (iread.eq.0) go to 10
1258 c write(*,*)'---> ',line
1264 iread=iread+ilen(var)+1
1265 read (line(iread:),*,err=12,end=12) val
1266 c write(*,*)'OK: ',var,' = ',val
1272 #elif (defined CRAY)
1273 integer lennam,lenval,ierror
1275 c getenv using a POSIX call, useful on the T3D
1276 c Sept 1996, comment out error check on advice of H. Pritchard
1279 if(lennam.le.0) stop '--error calling getenv--'
1280 call pxfgetenv(var,lennam,val,lenval,ierror)
1281 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1283 call getenv(var,val)