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,vbldpdum,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,vbldpdum,akp,rjunk,mp,ip,pstok
77 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
78 & j=1,nbondterm(i)),msc(i),isc(i),restok(i)
83 dsc_inv(i)=1.0D0/dsc(i)
88 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
89 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',
91 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
93 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
94 & vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
96 write (iout,'(13x,3f10.5)')
97 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
101 C reading lipid parameters
102 read(iliptranpar,*) pepliptran
104 read(iliptranpar,*) liptranene(i)
109 C Read the parameters of the probability distribution/energy expression
110 C of the virtual-bond valence angles theta
113 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),
114 & (bthet(j,i,1,1),j=1,2)
115 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
116 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
117 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
121 athet(1,i,1,-1)=athet(1,i,1,1)
122 athet(2,i,1,-1)=athet(2,i,1,1)
123 bthet(1,i,1,-1)=-bthet(1,i,1,1)
124 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)
132 athet(1,i,-1,-1)=athet(1,-i,1,1)
133 athet(2,i,-1,-1)=-athet(2,-i,1,1)
134 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
135 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
136 athet(1,i,-1,1)=athet(1,-i,1,1)
137 athet(2,i,-1,1)=-athet(2,-i,1,1)
138 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
139 bthet(2,i,-1,1)=bthet(2,-i,1,1)
140 athet(1,i,1,-1)=-athet(1,-i,1,1)
141 athet(2,i,1,-1)=athet(2,-i,1,1)
142 bthet(1,i,1,-1)=bthet(1,-i,1,1)
143 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
148 polthet(j,i)=polthet(j,-i)
151 gthet(j,i)=gthet(j,-i)
159 & 'Parameters of the virtual-bond valence angles:'
160 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
161 & ' ATHETA0 ',' A1 ',' A2 ',
164 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
165 & a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
167 write (iout,'(/a/9x,5a/79(1h-))')
168 & 'Parameters of the expression for sigma(theta_c):',
169 & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
170 & ' ALPH3 ',' SIGMA0C '
172 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
173 & (polthet(j,i),j=0,3),sigc0(i)
175 write (iout,'(/a/9x,5a/79(1h-))')
176 & 'Parameters of the second gaussian:',
177 & ' THETA0 ',' SIGMA0 ',' G1 ',
180 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
181 & sig0(i),(gthet(j,i),j=1,3)
185 & 'Parameters of the virtual-bond valence angles:'
186 write (iout,'(/a/9x,5a/79(1h-))')
187 & 'Coefficients of expansion',
188 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
189 & ' b1*10^1 ',' b2*10^1 '
191 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
192 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
193 & (10*bthet(j,i,1,1),j=1,2)
195 write (iout,'(/a/9x,5a/79(1h-))')
196 & 'Parameters of the expression for sigma(theta_c):',
197 & ' alpha0 ',' alph1 ',' alph2 ',
198 & ' alhp3 ',' sigma0c '
200 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
201 & (polthet(j,i),j=0,3),sigc0(i)
203 write (iout,'(/a/9x,5a/79(1h-))')
204 & 'Parameters of the second gaussian:',
205 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
208 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
209 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
215 C Read the parameters of Utheta determined from ab initio surfaces
216 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
218 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,
219 & ntheterm3,nsingle,ndouble
220 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
221 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
223 ithetyp(i)=-ithetyp(-i)
226 do i=-maxthetyp,maxthetyp
227 do j=-maxthetyp,maxthetyp
228 do k=-maxthetyp,maxthetyp
229 aa0thet(i,j,k,iblock)=0.0d0
231 aathet(l,i,j,k,iblock)=0.0d0
235 bbthet(m,l,i,j,k,iblock)=0.0d0
236 ccthet(m,l,i,j,k,iblock)=0.0d0
237 ddthet(m,l,i,j,k,iblock)=0.0d0
238 eethet(m,l,i,j,k,iblock)=0.0d0
244 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
245 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
253 c VAR:iblock means terminally blocking group 1=non-proline 2=proline
255 c VAR:ntethtyp is type of theta potentials type currently 0=glycine
256 c VAR:1=non-glicyne non-proline 2=proline
257 c VAR:negative values for D-aminoacid
259 do j=-nthetyp,nthetyp
260 do k=-nthetyp,nthetyp
261 read (ithep,'(6a)',end=111,err=111) res1
262 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
263 c VAR: aa0thet is variable describing the average value of Foureir
264 c VAR: expansion series
265 c VAR: aathet is foureir expansion in theta/2 angle for full formula
266 c VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
267 Condens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
268 read (ithep,*,end=111,err=111)
269 &(aathet(l,i,j,k,iblock),l=1,ntheterm)
270 read (ithep,*,end=111,err=111)
271 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
272 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
273 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
274 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),
276 read (ithep,*,end=111,err=111)
277 & (((ffthet(llll,lll,ll,i,j,k,iblock),
278 & ffthet(lll,llll,ll,i,j,k,iblock),
279 & ggthet(llll,lll,ll,i,j,k,iblock),
280 & ggthet(lll,llll,ll,i,j,k,iblock),
281 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
286 C For dummy ends assign glycine-type coefficients of theta-only terms; the
287 C coefficients of theta-and-gamma-dependent terms are zero.
288 C IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
289 C RECOMENTDED AFTER VERSION 3.3)
293 c aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
294 c aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
296 c aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
297 c aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
300 c aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
302 c aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
305 C AND COMMENT THE LOOPS BELOW
309 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
310 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
312 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
313 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
316 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
318 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
322 C Substitution for D aminoacids from symmetry.
325 do j=-nthetyp,nthetyp
326 do k=-nthetyp,nthetyp
327 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
329 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
333 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
334 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
335 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
336 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
342 ffthet(llll,lll,ll,i,j,k,iblock)=
343 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
344 ffthet(lll,llll,ll,i,j,k,iblock)=
345 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
346 ggthet(llll,lll,ll,i,j,k,iblock)=
347 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
348 ggthet(lll,llll,ll,i,j,k,iblock)=
349 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
358 C Control printout of the coefficients of virtual-bond-angle potentials
361 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
365 write (iout,'(//4a)')
366 & 'Type ',onelett(i),onelett(j),onelett(k)
367 write (iout,'(//a,10x,a)') " l","a[l]"
368 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
369 write (iout,'(i2,1pe15.5)')
370 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
372 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))')
373 & "b",l,"c",l,"d",l,"e",l
375 write (iout,'(i2,4(1pe15.5))') m,
376 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
377 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
381 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
382 & "f+",l,"f-",l,"g+",l,"g-",l
385 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
386 & ffthet(n,m,l,i,j,k,iblock),
387 & ffthet(m,n,l,i,j,k,iblock),
388 & ggthet(n,m,l,i,j,k,iblock),
389 & ggthet(m,n,l,i,j,k,iblock)
398 write (2,*) "Start reading THETA_PDB",ithep_pdb
401 read (ithep_pdb,*,err=111,end=111)
402 & a0thet(i),(athet(j,i,1,1),j=1,2),
403 & (bthet(j,i,1,1),j=1,2)
404 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
405 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
406 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
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)
414 athet(1,i,-1,1)=-athet(1,i,1,1)
415 athet(2,i,-1,1)=-athet(2,i,1,1)
416 bthet(1,i,-1,1)=bthet(1,i,1,1)
417 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)
429 athet(1,i,1,-1)=-athet(1,-i,1,1)
430 athet(2,i,1,-1)=athet(2,-i,1,1)
431 bthet(1,i,1,-1)=bthet(1,-i,1,1)
432 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
437 polthet(j,i)=polthet(j,-i)
440 gthet(j,i)=gthet(j,-i)
443 write (2,*) "End reading THETA_PDB"
449 C Read the parameters of the probability distribution/energy expression
450 C of the side chains.
453 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
457 dsc_inv(i)=1.0D0/dsc(i)
468 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
469 & ((blower(k,l,1),l=1,k),k=1,3)
470 censc(1,1,-i)=censc(1,1,i)
471 censc(2,1,-i)=censc(2,1,i)
472 censc(3,1,-i)=-censc(3,1,i)
474 read (irotam,*,end=112,err=112) bsc(j,i)
475 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
476 & ((blower(k,l,j),l=1,k),k=1,3)
477 censc(1,j,-i)=censc(1,j,i)
478 censc(2,j,-i)=censc(2,j,i)
479 censc(3,j,-i)=-censc(3,j,i)
480 C BSC is amplitude of Gaussian
487 akl=akl+blower(k,m,j)*blower(l,m,j)
491 if (((k.eq.3).and.(l.ne.3))
492 & .or.((l.eq.3).and.(k.ne.3))) then
493 gaussc(k,l,j,-i)=-akl
494 gaussc(l,k,j,-i)=-akl
506 write (iout,'(/a)') 'Parameters of side-chain local geometry'
511 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
512 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
513 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
514 & 'log h',(bsc(j,i),j=1,nlobi)
515 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
516 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
518 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
519 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
522 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
523 write (iout,'(a,f10.4,4(16x,f10.4))')
524 & 'Center ',(bsc(j,i),j=1,nlobi)
525 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
534 C Read scrot parameters for potentials determined from all-atom AM1 calculations
535 C added by Urszula Kozlowska 07/11/2007
538 read (irotam,*,end=112,err=112)
540 read (irotam,*,end=112,err=112)
543 read(irotam,*,end=112,err=112) sc_parmin(j,i)
548 C Read the parameters of the probability distribution/energy expression
549 C of the side chains.
551 write (2,*) "Start reading ROTAM_PDB"
553 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
557 dsc_inv(i)=1.0D0/dsc(i)
568 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
569 & ((blower(k,l,1),l=1,k),k=1,3)
571 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
572 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
573 & ((blower(k,l,j),l=1,k),k=1,3)
580 akl=akl+blower(k,m,j)*blower(l,m,j)
590 write (2,*) "End reading ROTAM_PDB"
596 C Read torsional parameters in old format
598 read (itorp,*,end=113,err=113) ntortyp,nterm_old
599 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
600 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
605 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
611 write (iout,'(/a/)') 'Torsional constants:'
614 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
615 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
621 C Read torsional parameters
623 read (itorp,*,end=113,err=113) ntortyp
624 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
627 itortyp(i)=-itortyp(-i)
629 write (iout,*) 'ntortyp',ntortyp
631 do j=-ntortyp+1,ntortyp-1
632 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
634 nterm(-i,-j,iblock)=nterm(i,j,iblock)
635 nlor(-i,-j,iblock)=nlor(i,j,iblock)
638 do k=1,nterm(i,j,iblock)
639 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
641 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
642 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
643 v0ij=v0ij+si*v1(k,i,j,iblock)
645 c write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
646 c write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
647 c &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
649 do k=1,nlor(i,j,iblock)
650 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
651 & vlor2(k,i,j),vlor3(k,i,j)
652 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
655 v0(-i,-j,iblock)=v0ij
661 write (iout,'(/a/)') 'Torsional constants:'
664 write (iout,*) 'ityp',i,' jtyp',j
665 write (iout,*) 'Fourier constants'
666 do k=1,nterm(i,j,iblock)
667 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
670 write (iout,*) 'Lorenz constants'
671 do k=1,nlor(i,j,iblock)
672 write (iout,'(3(1pe15.5))')
673 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
680 C 6/23/01 Read parameters for double torsionals
684 do j=-ntortyp+1,ntortyp-1
685 do k=-ntortyp+1,ntortyp-1
686 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
687 c write (iout,*) "OK onelett",
690 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
691 & .or. t3.ne.toronelet(k)) then
692 write (iout,*) "Error in double torsional parameter file",
695 call MPI_Finalize(Ierror)
697 stop "Error in double torsional parameter file"
699 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
700 & ntermd_2(i,j,k,iblock)
701 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
702 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
703 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
704 & ntermd_1(i,j,k,iblock))
705 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
706 & ntermd_1(i,j,k,iblock))
707 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
708 & ntermd_1(i,j,k,iblock))
709 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
710 & ntermd_1(i,j,k,iblock))
711 C Martix of D parameters for one dimesional foureir series
712 do l=1,ntermd_1(i,j,k,iblock)
713 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
714 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
715 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
716 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
717 c write(iout,*) "whcodze" ,
718 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
720 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
721 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
722 & v2s(m,l,i,j,k,iblock),
723 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
724 C Martix of D parameters for two dimesional fourier series
725 do l=1,ntermd_2(i,j,k,iblock)
727 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
728 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
729 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
730 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
739 write (iout,*) 'Constants for double torsionals'
742 do j=-ntortyp+1,ntortyp-1
743 do k=-ntortyp+1,ntortyp-1
744 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
745 & ' nsingle',ntermd_1(i,j,k,iblock),
746 & ' ndouble',ntermd_2(i,j,k,iblock)
748 write (iout,*) 'Single angles:'
749 do l=1,ntermd_1(i,j,k,iblock)
750 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
751 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
752 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
753 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
756 write (iout,*) 'Pairs of angles:'
757 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
758 do l=1,ntermd_2(i,j,k,iblock)
759 write (iout,'(i5,20f10.5)')
760 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
763 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
764 do l=1,ntermd_2(i,j,k,iblock)
765 write (iout,'(i5,20f10.5)')
766 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
767 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
776 C Read of Side-chain backbone correlation parameters
777 C Modified 11 May 2012 by Adasko
780 read (isccor,*,end=119,err=119) nsccortyp
782 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
784 isccortyp(i)=-isccortyp(-i)
786 iscprol=isccortyp(20)
787 c write (iout,*) 'ntortyp',ntortyp
789 cc maxinter is maximum interaction sites
793 read (isccor,*,end=119,err=119)
794 &nterm_sccor(i,j),nlor_sccor(i,j)
800 nterm_sccor(-i,j)=nterm_sccor(i,j)
801 nterm_sccor(-i,-j)=nterm_sccor(i,j)
802 nterm_sccor(i,-j)=nterm_sccor(i,j)
803 do k=1,nterm_sccor(i,j)
804 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
806 if (j.eq.iscprol) then
807 if (i.eq.isccortyp(10)) then
808 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
809 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
811 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
812 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
813 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
814 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
815 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
816 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
817 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
818 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
821 if (i.eq.isccortyp(10)) then
822 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
823 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
825 if (j.eq.isccortyp(10)) then
826 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
827 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
829 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
830 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
831 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
832 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
833 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
834 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
838 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
839 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
840 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
841 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
844 do k=1,nlor_sccor(i,j)
845 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
846 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
847 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
848 &(1+vlor3sccor(k,i,j)**2)
850 v0sccor(l,i,j)=v0ijsccor
851 v0sccor(l,-i,j)=v0ijsccor1
852 v0sccor(l,i,-j)=v0ijsccor2
853 v0sccor(l,-i,-j)=v0ijsccor3
859 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
860 c write (iout,*) 'ntortyp',ntortyp
862 cc maxinter is maximum interaction sites
866 read (isccor,*,end=119,err=119)
867 & nterm_sccor(i,j),nlor_sccor(i,j)
871 do k=1,nterm_sccor(i,j)
872 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
874 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
877 do k=1,nlor_sccor(i,j)
878 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
879 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
880 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
881 &(1+vlor3sccor(k,i,j)**2)
883 v0sccor(i,j,iblock)=v0ijsccor
891 write (iout,'(/a/)') 'Torsional constants:'
894 write (iout,*) 'ityp',i,' jtyp',j
895 write (iout,*) 'Fourier constants'
896 do k=1,nterm_sccor(i,j)
897 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
899 write (iout,*) 'Lorenz constants'
900 do k=1,nlor_sccor(i,j)
901 write (iout,'(3(1pe15.5))')
902 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
909 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
910 C interaction energy of the Gly, Ala, and Pro prototypes.
914 write (iout,*) "Coefficients of the cumulants"
916 read (ifourier,*) nloctyp
919 read (ifourier,*,end=115,err=115)
920 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
922 read (ifourier,*,end=115,err=115) (bnew1(ii,1,i),ii=1,3)
923 read (ifourier,*,end=115,err=115) (bnew2(ii,1,i),ii=1,3)
924 read (ifourier,*,end=115,err=115) (bnew1(ii,2,i),ii=1,1)
925 read (ifourier,*,end=115,err=115) (bnew2(ii,2,i),ii=1,1)
926 read (ifourier,*,end=115,err=115) (eenew(ii,i),ii=1,1)
929 write (iout,*) 'Type',i
930 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
938 c B1tilde(1,i) = b(3)
939 c B1tilde(2,i) =-b(5)
940 c B1tilde(1,-i) =-b(3)
941 c B1tilde(2,-i) =b(5)
975 Ctilde(2,1,i)=-b(9,i)
977 Ctilde(1,1,-i)=b(7,i)
978 Ctilde(1,2,-i)=-b(9,i)
979 Ctilde(2,1,-i)=b(9,i)
980 Ctilde(2,2,-i)=b(7,i)
982 c Ctilde(1,1,i)=0.0d0
983 c Ctilde(1,2,i)=0.0d0
984 c Ctilde(2,1,i)=0.0d0
985 c Ctilde(2,2,i)=0.0d0
1000 Dtilde(2,1,i)=-b(8,i)
1001 Dtilde(2,2,i)=b(6,i)
1002 Dtilde(1,1,-i)=b(6,i)
1003 Dtilde(1,2,-i)=-b(8,i)
1004 Dtilde(2,1,-i)=b(8,i)
1005 Dtilde(2,2,-i)=b(6,i)
1007 c Dtilde(1,1,i)=0.0d0
1008 c Dtilde(1,2,i)=0.0d0
1009 c Dtilde(2,1,i)=0.0d0
1010 c Dtilde(2,2,i)=0.0d0
1011 EEold(1,1,i)= b(10,i)+b(11,i)
1012 EEold(2,2,i)=-b(10,i)+b(11,i)
1013 EEold(2,1,i)= b(12,i)-b(13,i)
1014 EEold(1,2,i)= b(12,i)+b(13,i)
1015 EEold(1,1,-i)= b(10,i)+b(11,i)
1016 EEold(2,2,-i)=-b(10,i)+b(11,i)
1017 EEold(2,1,-i)=-b(12,i)+b(13,i)
1018 EEold(1,2,-i)=-b(12,i)-b(13,i)
1019 write(iout,*) "TU DOCHODZE"
1025 c ee(2,1,i)=ee(1,2,i)
1030 write (iout,*) 'Type',i
1032 write(iout,*) B1(1,i),B1(2,i)
1034 write(iout,*) B2(1,i),B2(2,i)
1037 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1041 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1045 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1052 C Read electrostatic-interaction parameters
1056 write (iout,'(/a)') 'Electrostatic interaction constants:'
1057 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1058 & 'IT','JT','APP','BPP','AEL6','AEL3'
1060 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1061 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1062 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1063 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1068 app (i,j)=epp(i,j)*rri*rri
1069 bpp (i,j)=-2.0D0*epp(i,j)*rri
1070 ael6(i,j)=elpp6(i,j)*4.2D0**6
1071 ael3(i,j)=elpp3(i,j)*4.2D0**3
1073 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1074 & ael6(i,j),ael3(i,j)
1079 C Read side-chain interaction parameters.
1081 read (isidep,*,end=117,err=117) ipot,expon
1082 if (ipot.lt.1 .or. ipot.gt.5) then
1083 write (iout,'(2a)') 'Error while reading SC interaction',
1084 & 'potential file - unknown potential type.'
1086 call MPI_Finalize(Ierror)
1092 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1093 & ', exponents are ',expon,2*expon
1094 goto (10,20,30,30,40) ipot
1095 C----------------------- LJ potential ---------------------------------
1096 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1097 & (sigma0(i),i=1,ntyp)
1099 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1100 write (iout,'(a/)') 'The epsilon array:'
1101 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1102 write (iout,'(/a)') 'One-body parameters:'
1103 write (iout,'(a,4x,a)') 'residue','sigma'
1104 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1107 C----------------------- LJK potential --------------------------------
1108 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1109 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1111 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1112 write (iout,'(a/)') 'The epsilon array:'
1113 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1114 write (iout,'(/a)') 'One-body parameters:'
1115 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1116 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1120 C---------------------- GB or BP potential -----------------------------
1122 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1124 read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
1125 read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
1126 read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
1127 read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
1128 C now we start reading lipid
1130 read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
1133 epslip(i,j)=epslip(i,j)+0.05d0
1136 C For the GB potential convert sigma'**2 into chi'
1139 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1143 write (iout,'(/a/)') 'Parameters of the BP potential:'
1144 write (iout,'(a/)') 'The epsilon array:'
1145 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1146 write (iout,'(/a)') 'One-body parameters:'
1147 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1149 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1150 & chip(i),alp(i),i=1,ntyp)
1153 C--------------------- GBV potential -----------------------------------
1154 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1155 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1156 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1158 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1159 write (iout,'(a/)') 'The epsilon array:'
1160 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1161 write (iout,'(/a)') 'One-body parameters:'
1162 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1163 & 's||/s_|_^2',' chip ',' alph '
1164 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1165 & sigii(i),chip(i),alp(i),i=1,ntyp)
1169 C-----------------------------------------------------------------------
1170 C Calculate the "working" parameters of SC interactions.
1174 epslip(i,j)=epslip(j,i)
1179 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1180 sigma(j,i)=sigma(i,j)
1181 rs0(i,j)=dwa16*sigma(i,j)
1185 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1186 & 'Working parameters of the SC interactions:',
1187 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1192 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1201 sigeps=dsign(1.0D0,epsij)
1203 aa_aq(i,j)=epsij*rrij*rrij
1204 bb_aq(i,j)=-sigeps*epsij*rrij
1205 aa_aq(j,i)=aa_aq(i,j)
1206 bb_aq(j,i)=bb_aq(i,j)
1207 epsijlip=epslip(i,j)
1208 sigeps=dsign(1.0D0,epsijlip)
1209 epsijlip=dabs(epsijlip)
1210 aa_lip(i,j)=epsijlip*rrij*rrij
1211 bb_lip(i,j)=-sigeps*epsijlip*rrij
1212 aa_lip(j,i)=aa_lip(i,j)
1213 bb_lip(j,i)=bb_lip(i,j)
1215 sigt1sq=sigma0(i)**2
1216 sigt2sq=sigma0(j)**2
1219 ratsig1=sigt2sq/sigt1sq
1220 ratsig2=1.0D0/ratsig1
1221 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1222 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1223 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1227 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1228 sigmaii(i,j)=rsum_max
1229 sigmaii(j,i)=rsum_max
1231 c sigmaii(i,j)=r0(i,j)
1232 c sigmaii(j,i)=r0(i,j)
1234 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1235 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1236 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1237 augm(i,j)=epsij*r_augm**(2*expon)
1238 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1245 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1246 & restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1247 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1253 C Define the SC-p interaction constants (hard-coded; old style)
1256 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1258 c aad(i,1)=0.3D0*4.0D0**12
1259 C Following line for constants currently implemented
1260 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1261 aad(i,1)=1.5D0*4.0D0**12
1262 c aad(i,1)=0.17D0*5.6D0**12
1264 C "Soft" SC-p repulsion
1266 C Following line for constants currently implemented
1267 c aad(i,1)=0.3D0*4.0D0**6
1268 C "Hard" SC-p repulsion
1269 bad(i,1)=3.0D0*4.0D0**6
1270 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1279 C 8/9/01 Read the SC-p interaction constants from file
1282 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1285 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1286 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1287 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1288 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1292 write (iout,*) "Parameters of SC-p interactions:"
1294 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1295 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1301 C Define the constants of the disulfide bridge
1305 c Old arbitrary potential - commented out.
1310 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1311 c energy surface of diethyl disulfide.
1312 c A. Liwo and U. Kozlowska, 11/24/03
1328 C if(me.eq.king) then
1329 C write (iout,'(/a)') "Disulfide bridge parameters:"
1330 C write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1331 C write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1332 C write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1333 C write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1337 111 write (iout,*) "Error reading bending energy parameters."
1339 112 write (iout,*) "Error reading rotamer energy parameters."
1341 113 write (iout,*) "Error reading torsional energy parameters."
1343 114 write (iout,*) "Error reading double torsional energy parameters."
1346 & "Error reading cumulant (multibody energy) parameters."
1348 116 write (iout,*) "Error reading electrostatic energy parameters."
1350 1161 write (iout,*) "Error reading electrostatic energy parameters.Lip"
1352 117 write (iout,*) "Error reading side chain interaction parameters."
1354 118 write (iout,*) "Error reading SCp interaction parameters."
1356 119 write (iout,*) "Error reading SCCOR parameters"
1359 call MPI_Finalize(Ierror)
1366 subroutine getenv_loc(var, val)
1367 character(*) var, val
1370 character(2000) line
1373 open (196,file='env',status='old',readonly,shared)
1375 c write(*,*)'looking for ',var
1376 10 read(196,*,err=11,end=11)line
1377 iread=index(line,var)
1378 c write(*,*)iread,' ',var,' ',line
1379 if (iread.eq.0) go to 10
1380 c write(*,*)'---> ',line
1386 iread=iread+ilen(var)+1
1387 read (line(iread:),*,err=12,end=12) val
1388 c write(*,*)'OK: ',var,' = ',val
1394 #elif (defined CRAY)
1395 integer lennam,lenval,ierror
1397 c getenv using a POSIX call, useful on the T3D
1398 c Sept 1996, comment out error check on advice of H. Pritchard
1401 if(lennam.le.0) stop '--error calling getenv--'
1402 call pxfgetenv(var,lennam,val,lenval,ierror)
1403 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1405 call getenv(var,val)