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 write(iout,*) "I am here",ntyp1
216 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
217 C write(iout,*) "I am herew"
219 ithetyp(i)=-ithetyp(-i)
222 do i=-maxthetyp,maxthetyp
223 do j=-maxthetyp,maxthetyp
224 do k=-maxthetyp,maxthetyp
225 aa0thet(i,j,k,iblock)=0.0d0
227 aathet(l,i,j,k,iblock)=0.0d0
231 bbthet(m,l,i,j,k,iblock)=0.0d0
232 ccthet(m,l,i,j,k,iblock)=0.0d0
233 ddthet(m,l,i,j,k,iblock)=0.0d0
234 eethet(m,l,i,j,k,iblock)=0.0d0
240 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
241 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
249 c VAR:iblock means terminally blocking group 1=non-proline 2=proline
251 c VAR:ntethtyp is type of theta potentials type currently 0=glycine
252 c VAR:1=non-glicyne non-proline 2=proline
253 c VAR:negative values for D-aminoacid
256 do j=-nthetyp,nthetyp
257 do k=-nthetyp,nthetyp
258 read (ithep,'(6a)',end=111,err=111) res1
259 c VAR: aa0thet is variable describing the average value of Foureir
260 c VAR: expansion series
261 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
262 c VAR: aathet is foureir expansion in theta/2 angle for full formula
263 c VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
264 Condens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
265 read (ithep,*,end=111,err=111)
266 &(aathet(l,i,j,k,iblock),l=1,ntheterm)
267 read (ithep,*,end=111,err=111)
268 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
269 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
270 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
271 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),
273 read (ithep,*,end=111,err=111)
274 & (((ffthet(llll,lll,ll,i,j,k,iblock),
275 & ffthet(lll,llll,ll,i,j,k,iblock),
276 & ggthet(llll,lll,ll,i,j,k,iblock),
277 & ggthet(lll,llll,ll,i,j,k,iblock),
278 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
285 C For dummy ends assign glycine-type coefficients of theta-only terms; the
286 C coefficients of theta-and-gamma-dependent terms are zero.
287 C IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
288 C RECOMENTDED AFTER VERSION 3.3)
292 c aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
293 c aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
295 c aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
296 c aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
299 c aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
301 c aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
304 C AND COMMENT THE LOOPS BELOW
308 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
309 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
311 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
312 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
315 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
317 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
321 C Substitution for D aminoacids from symmetry.
324 do j=-nthetyp,nthetyp
325 do k=-nthetyp,nthetyp
326 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
328 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
332 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
333 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
334 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
335 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
341 ffthet(llll,lll,ll,i,j,k,iblock)=
342 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
343 ffthet(lll,llll,ll,i,j,k,iblock)=
344 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
345 ggthet(llll,lll,ll,i,j,k,iblock)=
346 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
347 ggthet(lll,llll,ll,i,j,k,iblock)=
348 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
357 C Control printout of the coefficients of virtual-bond-angle potentials
360 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
364 write (iout,'(//4a)')
365 & 'Type ',onelett(i),onelett(j),onelett(k)
366 write (iout,'(//a,10x,a)') " l","a[l]"
367 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
368 write (iout,'(i2,1pe15.5)')
369 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
371 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))')
372 & "b",l,"c",l,"d",l,"e",l
374 write (iout,'(i2,4(1pe15.5))') m,
375 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
376 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
380 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
381 & "f+",l,"f-",l,"g+",l,"g-",l
384 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
385 & ffthet(n,m,l,i,j,k,iblock),
386 & ffthet(m,n,l,i,j,k,iblock),
387 & ggthet(n,m,l,i,j,k,iblock),
388 & ggthet(m,n,l,i,j,k,iblock)
397 write (2,*) "Start reading THETA_PDB"
399 read (ithep_pdb,*,err=111,end=111)
400 & a0thet(i),(athet(j,i,1,1),j=1,2),
401 & (bthet(j,i,1,1),j=1,2)
402 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
403 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
404 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
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)
412 athet(1,i,-1,1)=-athet(1,i,1,1)
413 athet(2,i,-1,1)=-athet(2,i,1,1)
414 bthet(1,i,-1,1)=bthet(1,i,1,1)
415 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)
427 athet(1,i,1,-1)=-athet(1,-i,1,1)
428 athet(2,i,1,-1)=athet(2,-i,1,1)
429 bthet(1,i,1,-1)=bthet(1,-i,1,1)
430 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
435 polthet(j,i)=polthet(j,-i)
438 gthet(j,i)=gthet(j,-i)
441 write (2,*) "End reading THETA_PDB"
447 C Read the parameters of the probability distribution/energy expression
448 C of the side chains.
451 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
455 dsc_inv(i)=1.0D0/dsc(i)
466 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
467 & ((blower(k,l,1),l=1,k),k=1,3)
469 read (irotam,*,end=112,err=112) bsc(j,i)
470 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
471 & ((blower(k,l,j),l=1,k),k=1,3)
478 akl=akl+blower(k,m,j)*blower(l,m,j)
489 write (iout,'(/a)') 'Parameters of side-chain local geometry'
494 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
495 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
496 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
497 & 'log h',(bsc(j,i),j=1,nlobi)
498 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
499 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
501 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
502 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
505 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
506 write (iout,'(a,f10.4,4(16x,f10.4))')
507 & 'Center ',(bsc(j,i),j=1,nlobi)
508 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
517 C Read scrot parameters for potentials determined from all-atom AM1 calculations
518 C added by Urszula Kozlowska 07/11/2007
521 read (irotam,*,end=112,err=112)
523 read (irotam,*,end=112,err=112)
526 read(irotam,*,end=112,err=112) sc_parmin(j,i)
531 C Read the parameters of the probability distribution/energy expression
532 C of the side chains.
534 write (2,*) "Start reading ROTAM_PDB"
537 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
541 dsc_inv(i)=1.0D0/dsc(i)
552 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
553 & ((blower(k,l,1),l=1,k),k=1,3)
555 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
556 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
557 & ((blower(k,l,j),l=1,k),k=1,3)
564 akl=akl+blower(k,m,j)*blower(l,m,j)
574 write (2,*) "Ending reading ROTAM_PDB"
580 C Read torsional parameters in old format
582 read (itorp,*,end=113,err=113) ntortyp,nterm_old
583 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
584 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
589 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
595 write (iout,'(/a/)') 'Torsional constants:'
598 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
599 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
605 C Read torsional parameters
607 read (itorp,*,end=113,err=113) ntortyp
608 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
611 itortyp(i)=-itortyp(-i)
613 write (iout,*) 'ntortyp',ntortyp
615 do j=-ntortyp+1,ntortyp-1
616 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
618 nterm(-i,-j,iblock)=nterm(i,j,iblock)
619 nlor(-i,-j,iblock)=nlor(i,j,iblock)
622 do k=1,nterm(i,j,iblock)
623 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
625 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
626 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
627 v0ij=v0ij+si*v1(k,i,j,iblock)
629 c write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
630 c write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
631 c &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
633 do k=1,nlor(i,j,iblock)
634 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
635 & vlor2(k,i,j),vlor3(k,i,j)
636 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
639 v0(-i,-j,iblock)=v0ij
645 write (iout,'(/a/)') 'Torsional constants:'
648 write (iout,*) 'ityp',i,' jtyp',j
649 write (iout,*) 'Fourier constants'
650 do k=1,nterm(i,j,iblock)
651 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
654 write (iout,*) 'Lorenz constants'
655 do k=1,nlor(i,j,iblock)
656 write (iout,'(3(1pe15.5))')
657 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
663 C 6/23/01 Read parameters for double torsionals
670 do j=-ntortyp+1,ntortyp-1
671 do k=-ntortyp+1,ntortyp-1
672 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
673 write (iout,*) "OK onelett",
676 if (t1.ne.onelett(i) .or. t2.ne.onelett(j)
677 & .or. t3.ne.onelett(k)) then
678 write (iout,*) "Error in double torsional parameter file",
681 call MPI_Finalize(Ierror)
683 stop "Error in double torsional parameter file"
685 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
686 & ntermd_2(i,j,k,iblock)
687 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
688 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
689 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
690 & ntermd_1(i,j,k,iblock))
691 read (itordp,*,end=114,err=114)(v1s(1,l,i,j,k,iblock),l=1,
692 & ntermd_1(i,j,k,iblock))
693 read (itordp,*,end=114,err=114)(v1c(2,l,i,j,k,iblock),l=1,
694 & ntermd_1(i,j,k,iblock))
695 read (itordp,*,end=114,err=114)(v1s(2,l,i,j,k,iblock),l=1,
696 & ntermd_1(i,j,k,iblock))
697 C Martix of D parameters for one dimesional foureir series
698 do l=1,ntermd_1(i,j,k,iblock)
699 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
700 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
701 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
702 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
703 c write(iout,*) "whcodze" ,
704 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
706 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
707 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
708 & v2s(m,l,i,j,k,iblock),
709 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
710 C Martix of D parameters for two dimesional fourier series
711 do l=1,ntermd_2(i,j,k,iblock)
713 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
714 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
715 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
716 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
725 write (iout,*) 'Constants for double torsionals'
728 do j=-ntortyp+1,ntortyp-1
729 do k=-ntortyp+1,ntortyp-1
730 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
731 & ' nsingle',ntermd_1(i,j,k,iblock),
732 & ' ndouble',ntermd_2(i,j,k,iblock)
734 write (iout,*) 'Single angles:'
735 do l=1,ntermd_1(i,j,k,iblock)
736 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
737 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
738 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
739 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
742 write (iout,*) 'Pairs of angles:'
743 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
744 do l=1,ntermd_2(i,j,k,iblock)
745 write (iout,'(i5,20f10.5)')
746 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
749 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
750 do l=1,ntermd_2(i,j,k,iblock)
751 write (iout,'(i5,20f10.5)')
752 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
753 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
762 C Read of Side-chain backbone correlation parameters
763 C Modified 11 May 2012 by Adasko
766 read (isccor,*,end=1113,err=1113) nsccortyp
768 write (iout,*) "Tu wchodze"
769 read (isccor,*,end=1113,err=1113) (isccortyp(i),i=1,ntyp)
771 isccortyp(i)=-isccortyp(-i)
773 iscprol=isccortyp(20)
774 C write (iout,*) 'ntortyp',ntortyp
776 cc maxinter is maximum interaction sites
780 read (isccor,*,end=1113,err=1113) nterm_sccor(i,j),
788 nterm_sccor(-i,j)=nterm_sccor(i,j)
789 nterm_sccor(-i,-j)=nterm_sccor(i,j)
790 nterm_sccor(i,-j)=nterm_sccor(i,j)
791 do k=1,nterm_sccor(i,j)
792 read (isccor,*,end=1113,err=1113) kk,v1sccor(k,l,i,j)
794 c write(iout,*) "k=",kk
795 if (j.eq.iscprol) then
796 if (i.eq.isccortyp(10)) then
797 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
798 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
800 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
801 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
802 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
803 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
804 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
805 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
806 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
807 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
810 if (i.eq.isccortyp(10)) then
811 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
812 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
814 if (j.eq.isccortyp(10)) then
815 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
816 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
818 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
819 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
820 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
821 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
822 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
823 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
827 C read (isccor,*,end=1113,err=1113) kk,v1sccor(k,l,i,j)
828 C & ,v2sccor(k,l,i,j)
829 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
830 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
831 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
832 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
835 do k=1,nlor_sccor(i,j)
836 read (isccor,*,end=1113,err=1113) kk,vlor1sccor(k,i,j),
837 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
838 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
839 &(1+vlor3sccor(k,i,j)**2)
841 v0sccor(l,i,j)=v0ijsccor
842 v0sccor(l,-i,j)=v0ijsccor1
843 v0sccor(l,i,-j)=v0ijsccor2
844 v0sccor(l,-i,-j)=v0ijsccor3
850 write(iout,*) "a tu nie wchodze"
851 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
852 c write (iout,*) 'ntortyp',ntortyp
854 cc maxinter is maximum interaction sites
858 read (isccor,*,end=119,err=119)
859 & nterm_sccor(i,j),nlor_sccor(i,j)
863 do k=1,nterm_sccor(i,j)
864 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
866 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
869 do k=1,nlor_sccor(i,j)
870 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
871 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
872 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
873 &(1+vlor3sccor(k,i,j)**2)
875 v0sccor(i,j,iblock)=v0ijsccor
883 write (iout,'(/a/)') 'Torsional constants:'
886 write (iout,*) 'ityp',i,' jtyp',j
887 write (iout,*) 'Fourier constants'
888 do k=1,nterm_sccor(i,j)
889 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
891 write (iout,*) 'Lorenz constants'
892 do k=1,nlor_sccor(i,j)
893 write (iout,'(3(1pe15.5))')
894 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
901 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
902 C interaction energy of the Gly, Ala, and Pro prototypes.
906 write (iout,*) "Coefficients of the cumulants"
908 read (ifourier,*) nloctyp
910 read (ifourier,*,end=115,err=115)
911 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
913 write (iout,*) 'Type',i
914 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
922 C B1tilde(1,-i) =-b(3)
923 C B1tilde(2,-i) =b(5)
942 c Ctilde(1,1,i)=0.0d0
943 c Ctilde(1,2,i)=0.0d0
944 c Ctilde(2,1,i)=0.0d0
945 c Ctilde(2,2,i)=0.0d0
958 c Dtilde(1,1,i)=0.0d0
959 c Dtilde(1,2,i)=0.0d0
960 c Dtilde(2,1,i)=0.0d0
961 c Dtilde(2,2,i)=0.0d0
962 EE(1,1,i)= b(10)+b(11)
963 EE(2,2,i)=-b(10)+b(11)
964 EE(2,1,i)= b(12)-b(13)
965 EE(1,2,i)= b(12)+b(13)
970 c ee(2,1,i)=ee(1,2,i)
974 write (iout,*) 'Type',i
976 write(iout,*) B1(1,i),B1(2,i)
978 write(iout,*) B2(1,i),B2(2,i)
981 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
985 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
989 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
994 C Read electrostatic-interaction parameters
998 write (iout,'(/a)') 'Electrostatic interaction constants:'
999 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1000 & 'IT','JT','APP','BPP','AEL6','AEL3'
1002 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1003 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1004 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1005 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1010 app (i,j)=epp(i,j)*rri*rri
1011 bpp (i,j)=-2.0D0*epp(i,j)*rri
1012 ael6(i,j)=elpp6(i,j)*4.2D0**6
1013 ael3(i,j)=elpp3(i,j)*4.2D0**3
1014 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1015 & ael6(i,j),ael3(i,j)
1019 C Read side-chain interaction parameters.
1021 read (isidep,*,end=117,err=117) ipot,expon
1022 if (ipot.lt.1 .or. ipot.gt.5) then
1023 write (iout,'(2a)') 'Error while reading SC interaction',
1024 & 'potential file - unknown potential type.'
1026 call MPI_Finalize(Ierror)
1032 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1033 & ', exponents are ',expon,2*expon
1034 goto (10,20,30,30,40) ipot
1035 C----------------------- LJ potential ---------------------------------
1036 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1037 & (sigma0(i),i=1,ntyp)
1039 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1040 write (iout,'(a/)') 'The epsilon array:'
1041 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1042 write (iout,'(/a)') 'One-body parameters:'
1043 write (iout,'(a,4x,a)') 'residue','sigma'
1044 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1047 C----------------------- LJK potential --------------------------------
1048 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1049 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1051 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1052 write (iout,'(a/)') 'The epsilon array:'
1053 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1054 write (iout,'(/a)') 'One-body parameters:'
1055 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1056 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1060 C---------------------- GB or BP potential -----------------------------
1062 read (isidep,*,end=116,err=116)(eps(i,j),j=i,ntyp)
1064 read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
1065 read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
1066 read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
1067 read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
1069 c 30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1070 c & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
1071 c & (alp(i),i=1,ntyp)
1072 C For the GB potential convert sigma'**2 into chi'
1075 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1079 write (iout,'(/a/)') 'Parameters of the BP potential:'
1080 write (iout,'(a/)') 'The epsilon array:'
1081 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1082 write (iout,'(/a)') 'One-body parameters:'
1083 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1085 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1086 & chip(i),alp(i),i=1,ntyp)
1089 C--------------------- GBV potential -----------------------------------
1090 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1091 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1092 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1094 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1095 write (iout,'(a/)') 'The epsilon array:'
1096 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1097 write (iout,'(/a)') 'One-body parameters:'
1098 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1099 & 's||/s_|_^2',' chip ',' alph '
1100 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1101 & sigii(i),chip(i),alp(i),i=1,ntyp)
1105 C-----------------------------------------------------------------------
1106 C Calculate the "working" parameters of SC interactions.
1114 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1115 sigma(j,i)=sigma(i,j)
1116 rs0(i,j)=dwa16*sigma(i,j)
1120 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1121 & 'Working parameters of the SC interactions:',
1122 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1127 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1136 sigeps=dsign(1.0D0,epsij)
1138 aa(i,j)=epsij*rrij*rrij
1139 bb(i,j)=-sigeps*epsij*rrij
1143 sigt1sq=sigma0(i)**2
1144 sigt2sq=sigma0(j)**2
1147 ratsig1=sigt2sq/sigt1sq
1148 ratsig2=1.0D0/ratsig1
1149 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1150 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1151 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1155 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1156 sigmaii(i,j)=rsum_max
1157 sigmaii(j,i)=rsum_max
1159 c sigmaii(i,j)=r0(i,j)
1160 c sigmaii(j,i)=r0(i,j)
1162 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1163 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1164 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1165 augm(i,j)=epsij*r_augm**(2*expon)
1166 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1173 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1174 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1175 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1181 C Define the SC-p interaction constants (hard-coded; old style)
1184 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1186 c aad(i,1)=0.3D0*4.0D0**12
1187 C Following line for constants currently implemented
1188 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1189 aad(i,1)=1.5D0*4.0D0**12
1190 c aad(i,1)=0.17D0*5.6D0**12
1192 C "Soft" SC-p repulsion
1194 C Following line for constants currently implemented
1195 c aad(i,1)=0.3D0*4.0D0**6
1196 C "Hard" SC-p repulsion
1197 bad(i,1)=3.0D0*4.0D0**6
1198 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1207 C 8/9/01 Read the SC-p interaction constants from file
1210 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1213 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1214 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1215 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1216 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1220 write (iout,*) "Parameters of SC-p interactions:"
1222 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1223 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1228 C Define the constants of the disulfide bridge
1232 c Old arbitrary potential - commented out.
1237 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1238 c energy surface of diethyl disulfide.
1239 c A. Liwo and U. Kozlowska, 11/24/03
1256 write (iout,'(/a)') "Disulfide bridge parameters:"
1257 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1258 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1259 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1260 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1264 111 write (iout,*) "Error reading bending energy parameters."
1266 112 write (iout,*) "Error reading rotamer energy parameters."
1268 113 write (iout,*) "Error reading torsional energy parameters."
1271 & "Error reading side-chain torsional energy parameters."
1273 114 write (iout,*) "Error reading double torsional energy parameters."
1276 & "Error reading cumulant (multibody energy) parameters."
1278 116 write (iout,*) "Error reading electrostatic energy parameters."
1280 117 write (iout,*) "Error reading side chain interaction parameters."
1282 118 write (iout,*) "Error reading SCp interaction parameters."
1284 119 write (iout,*) "Error reading SCCOR parameters"
1287 call MPI_Finalize(Ierror)
1294 subroutine getenv_loc(var, val)
1295 character(*) var, val
1298 character(2000) line
1301 open (196,file='env',status='old',readonly,shared)
1303 c write(*,*)'looking for ',var
1304 10 read(196,*,err=11,end=11)line
1305 iread=index(line,var)
1306 c write(*,*)iread,' ',var,' ',line
1307 if (iread.eq.0) go to 10
1308 c write(*,*)'---> ',line
1314 iread=iread+ilen(var)+1
1315 read (line(iread:),*,err=12,end=12) val
1316 c write(*,*)'OK: ',var,' = ',val
1322 #elif (defined CRAY)
1323 integer lennam,lenval,ierror
1325 c getenv using a POSIX call, useful on the T3D
1326 c Sept 1996, comment out error check on advice of H. Pritchard
1329 if(lennam.le.0) stop '--error calling getenv--'
1330 call pxfgetenv(var,lennam,val,lenval,ierror)
1331 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1333 call getenv(var,val)