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"/
32 dimension blower(3,3,maxlob)
34 character*3 lancuch,ucase
36 C For printing parameters after they are read set the following in the UNRES
39 C setenv PRINT_PARM YES
41 C To print parameters in LaTeX format rather than as ASCII tables:
45 call getenv_loc("PRINT_PARM",lancuch)
46 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
47 call getenv_loc("LATEX",lancuch)
48 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
50 dwa16=2.0d0**(1.0d0/6.0d0)
52 C Assign virtual-bond length
57 c Read the virtual-bond parameters, masses, and moments of inertia
58 c and Stokes' radii of the peptide group and side chains
61 read (ibond,*) vbldp0,akp,mp,ip,pstok
64 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
69 dsc_inv(i)=1.0D0/dsc(i)
73 read (ibond,*) junk,vbldp0,akp,rjunk,mp,ip,pstok
75 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
76 & j=1,nbondterm(i)),msc(i),isc(i),restok(i)
81 dsc_inv(i)=1.0D0/dsc(i)
86 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
87 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',
89 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
91 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
92 & vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
94 write (iout,'(13x,3f10.5)')
95 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
101 C Read the parameters of the probability distribution/energy expression
102 C of the virtual-bond valence angles theta
105 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i),j=1,2),
107 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
108 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
109 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
116 & 'Parameters of the virtual-bond valence angles:'
117 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
118 & ' ATHETA0 ',' A1 ',' A2 ',
121 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
122 & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
124 write (iout,'(/a/9x,5a/79(1h-))')
125 & 'Parameters of the expression for sigma(theta_c):',
126 & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
127 & ' ALPH3 ',' SIGMA0C '
129 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
130 & (polthet(j,i),j=0,3),sigc0(i)
132 write (iout,'(/a/9x,5a/79(1h-))')
133 & 'Parameters of the second gaussian:',
134 & ' THETA0 ',' SIGMA0 ',' G1 ',
137 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
138 & sig0(i),(gthet(j,i),j=1,3)
142 & 'Parameters of the virtual-bond valence angles:'
143 write (iout,'(/a/9x,5a/79(1h-))')
144 & 'Coefficients of expansion',
145 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
146 & ' b1*10^1 ',' b2*10^1 '
148 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
149 & a0thet(i),(100*athet(j,i),j=1,2),(10*bthet(j,i),j=1,2)
151 write (iout,'(/a/9x,5a/79(1h-))')
152 & 'Parameters of the expression for sigma(theta_c):',
153 & ' alpha0 ',' alph1 ',' alph2 ',
154 & ' alhp3 ',' sigma0c '
156 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
157 & (polthet(j,i),j=0,3),sigc0(i)
159 write (iout,'(/a/9x,5a/79(1h-))')
160 & 'Parameters of the second gaussian:',
161 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
164 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
165 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
171 C Read the parameters of Utheta determined from ab initio surfaces
172 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
174 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,
175 & ntheterm3,nsingle,ndouble
176 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
177 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
179 ithetyp(i)=-ithetyp(-i)
182 do i=-maxthetyp,maxthetyp
183 do j=-maxthetyp,maxthetyp
184 do k=-maxthetyp,maxthetyp
185 aa0thet(i,j,k,iblock)=0.0d0
187 aathet(l,i,j,k,iblock)=0.0d0
191 bbthet(m,l,i,j,k,iblock)=0.0d0
192 ccthet(m,l,i,j,k,iblock)=0.0d0
193 ddthet(m,l,i,j,k,iblock)=0.0d0
194 eethet(m,l,i,j,k,iblock)=0.0d0
200 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
201 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
209 c VAR:iblock means terminally blocking group 1=non-proline 2=proline
211 c VAR:ntethtyp is type of theta potentials type currently 0=glycine
212 c VAR:1=non-glicyne non-proline 2=proline
213 c VAR:negative values for D-aminoacid
215 do j=-nthetyp,nthetyp
216 do k=-nthetyp,nthetyp
217 read (ithep,'(6a)',end=111,err=111) res1
218 c VAR: aa0thet is variable describing the average value of Foureir
219 c VAR: expansion series
220 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
221 c VAR: aathet is foureir expansion in theta/2 angle for full formula
222 c VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
223 Condens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
224 read (ithep,*,end=111,err=111)
225 &(aathet(l,i,j,k,iblock),l=1,ntheterm)
226 read (ithep,*,end=111,err=111)
227 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
228 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
229 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
230 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),
232 read (ithep,*,end=111,err=111)
233 & (((ffthet(llll,lll,ll,i,j,k,iblock),
234 & ffthet(lll,llll,ll,i,j,k,iblock),
235 & ggthet(llll,lll,ll,i,j,k,iblock),
236 & ggthet(lll,llll,ll,i,j,k,iblock),
237 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
244 C For dummy ends assign glycine-type coefficients of theta-only terms; the
245 C coefficients of theta-and-gamma-dependent terms are zero.
246 C IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
247 C RECOMENTDED AFTER VERSION 3.3)
251 c aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
252 c aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
254 c aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
255 c aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
258 c aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
260 c aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
263 C AND COMMENT THE LOOPS BELOW
267 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
268 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
270 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
271 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
274 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
276 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
280 C Substitution for D aminoacids from symmetry.
283 do j=-nthetyp,nthetyp
284 do k=-nthetyp,nthetyp
285 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
287 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
291 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
292 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
293 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
294 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
300 ffthet(llll,lll,ll,i,j,k,iblock)=
301 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
302 ffthet(lll,llll,ll,i,j,k,iblock)=
303 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
304 ggthet(llll,lll,ll,i,j,k,iblock)=
305 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
306 ggthet(lll,llll,ll,i,j,k,iblock)=
307 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
316 C Control printout of the coefficients of virtual-bond-angle potentials
319 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
323 write (iout,'(//4a)')
324 & 'Type ',onelett(i),onelett(j),onelett(k)
325 write (iout,'(//a,10x,a)') " l","a[l]"
326 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
327 write (iout,'(i2,1pe15.5)')
328 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
330 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))')
331 & "b",l,"c",l,"d",l,"e",l
333 write (iout,'(i2,4(1pe15.5))') m,
334 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
335 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
339 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
340 & "f+",l,"f-",l,"g+",l,"g-",l
343 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
344 & ffthet(n,m,l,i,j,k,iblock),
345 & ffthet(m,n,l,i,j,k,iblock),
346 & ggthet(n,m,l,i,j,k,iblock),
347 & ggthet(m,n,l,i,j,k,iblock)
356 write (2,*) "Start reading THETA_PDB"
358 read (ithep_pdb,*,err=111,end=111)
359 & a0thet(i),(athet(j,i,1,1),j=1,2),
360 & (bthet(j,i,1,1),j=1,2)
361 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
362 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
363 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
367 athet(1,i,1,-1)=athet(1,i,1,1)
368 athet(2,i,1,-1)=athet(2,i,1,1)
369 bthet(1,i,1,-1)=-bthet(1,i,1,1)
370 bthet(2,i,1,-1)=-bthet(2,i,1,1)
371 athet(1,i,-1,1)=-athet(1,i,1,1)
372 athet(2,i,-1,1)=-athet(2,i,1,1)
373 bthet(1,i,-1,1)=bthet(1,i,1,1)
374 bthet(2,i,-1,1)=bthet(2,i,1,1)
378 athet(1,i,-1,-1)=athet(1,-i,1,1)
379 athet(2,i,-1,-1)=-athet(2,-i,1,1)
380 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
381 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
382 athet(1,i,-1,1)=athet(1,-i,1,1)
383 athet(2,i,-1,1)=-athet(2,-i,1,1)
384 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
385 bthet(2,i,-1,1)=bthet(2,-i,1,1)
386 athet(1,i,1,-1)=-athet(1,-i,1,1)
387 athet(2,i,1,-1)=athet(2,-i,1,1)
388 bthet(1,i,1,-1)=bthet(1,-i,1,1)
389 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
394 polthet(j,i)=polthet(j,-i)
397 gthet(j,i)=gthet(j,-i)
400 write (2,*) "End reading THETA_PDB"
406 C Read the parameters of the probability distribution/energy expression
407 C of the side chains.
410 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
414 dsc_inv(i)=1.0D0/dsc(i)
425 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
426 & ((blower(k,l,1),l=1,k),k=1,3)
428 read (irotam,*,end=112,err=112) bsc(j,i)
429 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
430 & ((blower(k,l,j),l=1,k),k=1,3)
437 akl=akl+blower(k,m,j)*blower(l,m,j)
448 write (iout,'(/a)') 'Parameters of side-chain local geometry'
453 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
454 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
455 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
456 & 'log h',(bsc(j,i),j=1,nlobi)
457 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
458 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
460 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
461 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
464 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
465 write (iout,'(a,f10.4,4(16x,f10.4))')
466 & 'Center ',(bsc(j,i),j=1,nlobi)
467 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
476 C Read scrot parameters for potentials determined from all-atom AM1 calculations
477 C added by Urszula Kozlowska 07/11/2007
480 read (irotam,*,end=112,err=112)
482 read (irotam,*,end=112,err=112)
485 read(irotam,*,end=112,err=112) sc_parmin(j,i)
490 C Read the parameters of the probability distribution/energy expression
491 C of the side chains.
493 write (2,*) "Start reading ROTAM_PDB"
496 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
500 dsc_inv(i)=1.0D0/dsc(i)
511 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
512 & ((blower(k,l,1),l=1,k),k=1,3)
514 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
515 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
516 & ((blower(k,l,j),l=1,k),k=1,3)
523 akl=akl+blower(k,m,j)*blower(l,m,j)
533 write (2,*) "Ending reading ROTAM_PDB"
539 C Read torsional parameters in old format
541 read (itorp,*,end=113,err=113) ntortyp,nterm_old
542 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
543 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
548 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
554 write (iout,'(/a/)') 'Torsional constants:'
557 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
558 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
564 C Read torsional parameters
566 read (itorp,*,end=113,err=113) ntortyp
567 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
568 c write (iout,*) 'ntortyp',ntortyp
571 read (itorp,*,end=113,err=113) nterm(i,j),nlor(i,j)
575 read (itorp,*,end=113,err=113) kk,v1(k,i,j),v2(k,i,j)
576 v0ij=v0ij+si*v1(k,i,j)
580 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
581 & vlor2(k,i,j),vlor3(k,i,j)
582 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
589 write (iout,'(/a/)') 'Torsional constants:'
592 write (iout,*) 'ityp',i,' jtyp',j
593 write (iout,*) 'Fourier constants'
595 write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
597 write (iout,*) 'Lorenz constants'
599 write (iout,'(3(1pe15.5))')
600 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
606 C 6/23/01 Read parameters for double torsionals
611 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
613 c write (iout,*) "OK onelett",
616 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
617 & .or. t3.ne.toronelet(k)) then
619 if (t1.ne.onelett(i) .or. t2.ne.onelett(j)
620 & .or. t3.ne.onelett(k)) then
622 write (iout,*) "Error in double torsional parameter file",
625 call MPI_Finalize(Ierror)
627 stop "Error in double torsional parameter file"
629 read (itordp,*,end=114,err=114) ntermd_1(i,j,k),
631 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k),l=1,
633 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k),l=1,
635 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k),l=1,
637 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k),l=1,
639 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k),
640 & v2c(m,l,i,j,k),v2s(l,m,i,j,k),v2s(m,l,i,j,k),
641 & m=1,l-1),l=1,ntermd_2(i,j,k))
647 write (iout,*) 'Constants for double torsionals'
651 do j=-ntortyp+1,ntortyp-1
652 do k=-ntortyp+1,ntortyp-1
658 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
659 & ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k)
661 write (iout,*) 'Single angles:'
662 do l=1,ntermd_1(i,j,k)
663 write (iout,'(i5,2f10.5,5x,2f10.5)') l,
664 & v1c(1,l,i,j,k),v1s(1,l,i,j,k),
665 & v1c(2,l,i,j,k),v1s(2,l,i,j,k)
668 write (iout,*) 'Pairs of angles:'
669 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
670 do l=1,ntermd_2(i,j,k)
671 write (iout,'(i5,20f10.5)')
672 & l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k))
675 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
676 do l=1,ntermd_2(i,j,k)
677 write (iout,'(i5,20f10.5)')
678 & l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k))
686 C Read of Side-chain backbone correlation parameters
687 C Modified 11 May 2012 by Adasko
691 read (isccor,*,end=119,err=119) nsccortyp
693 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
695 isccortyp(i)=-isccortyp(-i)
697 iscprol=isccortyp(20)
699 read (isccor,*,end=1113,err=1113) nsccortyp
700 read (isccor,*,end=1113,err=1113) (isccortyp(i),i=1,ntyp)
702 c write (iout,*) 'ntortyp',ntortyp
704 cc maxinter is maximum interaction sites
709 read (isccor,*,end=119,err=119) nterm_sccor(i,j),nlor_sccor(i,j)
711 read (isccor,*,end=1113,err=1113) nterm_sccor(i,j),
719 nterm_sccor(-i,j)=nterm_sccor(i,j)
720 nterm_sccor(-i,-j)=nterm_sccor(i,j)
721 nterm_sccor(i,-j)=nterm_sccor(i,j)
722 do k=1,nterm_sccor(i,j)
724 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
726 if (j.eq.iscprol) then
727 if (i.eq.isccortyp(10)) then
728 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
729 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
731 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
732 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
733 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
734 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
735 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
736 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
737 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
738 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
741 if (i.eq.isccortyp(10)) then
742 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
743 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
745 if (j.eq.isccortyp(10)) then
746 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
747 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
749 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
750 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
751 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
752 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
753 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
754 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
759 read (isccor,*,end=1113,err=1113) kk,v1sccor(k,l,i,j)
762 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
763 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
764 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
765 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
768 do k=1,nlor_sccor(i,j)
770 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
772 read (isccor,*,end=1113,err=1113) kk,vlor1sccor(k,i,j),
774 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
775 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
776 &(1+vlor3sccor(k,i,j)**2)
778 v0sccor(l,i,j)=v0ijsccor
779 v0sccor(l,-i,j)=v0ijsccor1
780 v0sccor(l,i,-j)=v0ijsccor2
781 v0sccor(l,-i,-j)=v0ijsccor3
787 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
788 c write (iout,*) 'ntortyp',ntortyp
790 cc maxinter is maximum interaction sites
794 read (isccor,*,end=119,err=119)
795 & nterm_sccor(i,j),nlor_sccor(i,j)
799 do k=1,nterm_sccor(i,j)
800 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
802 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
805 do k=1,nlor_sccor(i,j)
806 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
807 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
808 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
809 &(1+vlor3sccor(k,i,j)**2)
811 v0sccor(i,j,iblock)=v0ijsccor
819 write (iout,'(/a/)') 'Torsional constants:'
822 write (iout,*) 'ityp',i,' jtyp',j
823 write (iout,*) 'Fourier constants'
824 do k=1,nterm_sccor(i,j)
825 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
827 write (iout,*) 'Lorenz constants'
828 do k=1,nlor_sccor(i,j)
829 write (iout,'(3(1pe15.5))')
830 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
837 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
838 C interaction energy of the Gly, Ala, and Pro prototypes.
842 write (iout,*) "Coefficients of the cumulants"
844 read (ifourier,*) nloctyp
846 read (ifourier,*,end=115,err=115)
847 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
849 write (iout,*) 'Type',i
850 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
882 c Ctilde(1,1,i)=0.0d0
883 c Ctilde(1,2,i)=0.0d0
884 c Ctilde(2,1,i)=0.0d0
885 c Ctilde(2,2,i)=0.0d0
898 c Dtilde(1,1,i)=0.0d0
899 c Dtilde(1,2,i)=0.0d0
900 c Dtilde(2,1,i)=0.0d0
901 c Dtilde(2,2,i)=0.0d0
902 EE(1,1,i)= b(10)+b(11)
903 EE(2,2,i)=-b(10)+b(11)
904 EE(2,1,i)= b(12)-b(13)
905 EE(1,2,i)= b(12)+b(13)
910 c ee(2,1,i)=ee(1,2,i)
914 write (iout,*) 'Type',i
916 write(iout,*) B1(1,i),B1(2,i)
918 write(iout,*) B2(1,i),B2(2,i)
921 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
925 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
929 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
934 C Read electrostatic-interaction parameters
938 write (iout,'(/a)') 'Electrostatic interaction constants:'
939 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
940 & 'IT','JT','APP','BPP','AEL6','AEL3'
942 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
943 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
944 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
945 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
950 app (i,j)=epp(i,j)*rri*rri
951 bpp (i,j)=-2.0D0*epp(i,j)*rri
952 ael6(i,j)=elpp6(i,j)*4.2D0**6
953 ael3(i,j)=elpp3(i,j)*4.2D0**3
954 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
955 & ael6(i,j),ael3(i,j)
959 C Read side-chain interaction parameters.
961 read (isidep,*,end=117,err=117) ipot,expon
962 if (ipot.lt.1 .or. ipot.gt.5) then
963 write (iout,'(2a)') 'Error while reading SC interaction',
964 & 'potential file - unknown potential type.'
966 call MPI_Finalize(Ierror)
972 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
973 & ', exponents are ',expon,2*expon
974 goto (10,20,30,30,40) ipot
975 C----------------------- LJ potential ---------------------------------
976 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
977 & (sigma0(i),i=1,ntyp)
979 write (iout,'(/a/)') 'Parameters of the LJ potential:'
980 write (iout,'(a/)') 'The epsilon array:'
981 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
982 write (iout,'(/a)') 'One-body parameters:'
983 write (iout,'(a,4x,a)') 'residue','sigma'
984 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
987 C----------------------- LJK potential --------------------------------
988 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
989 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
991 write (iout,'(/a/)') 'Parameters of the LJK potential:'
992 write (iout,'(a/)') 'The epsilon array:'
993 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
994 write (iout,'(/a)') 'One-body parameters:'
995 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
996 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1000 C---------------------- GB or BP potential -----------------------------
1002 read (isidep,*,end=116,err=116)(eps(i,j),j=i,ntyp)
1004 read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
1005 read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
1006 read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
1007 read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
1009 c 30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1010 c & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
1011 c & (alp(i),i=1,ntyp)
1012 C For the GB potential convert sigma'**2 into chi'
1015 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1019 write (iout,'(/a/)') 'Parameters of the BP potential:'
1020 write (iout,'(a/)') 'The epsilon array:'
1021 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1022 write (iout,'(/a)') 'One-body parameters:'
1023 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1025 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1026 & chip(i),alp(i),i=1,ntyp)
1029 C--------------------- GBV potential -----------------------------------
1030 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1031 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1032 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1034 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1035 write (iout,'(a/)') 'The epsilon array:'
1036 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1037 write (iout,'(/a)') 'One-body parameters:'
1038 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1039 & 's||/s_|_^2',' chip ',' alph '
1040 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1041 & sigii(i),chip(i),alp(i),i=1,ntyp)
1045 C-----------------------------------------------------------------------
1046 C Calculate the "working" parameters of SC interactions.
1054 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1055 sigma(j,i)=sigma(i,j)
1056 rs0(i,j)=dwa16*sigma(i,j)
1060 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1061 & 'Working parameters of the SC interactions:',
1062 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1067 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1076 sigeps=dsign(1.0D0,epsij)
1078 aa(i,j)=epsij*rrij*rrij
1079 bb(i,j)=-sigeps*epsij*rrij
1083 sigt1sq=sigma0(i)**2
1084 sigt2sq=sigma0(j)**2
1087 ratsig1=sigt2sq/sigt1sq
1088 ratsig2=1.0D0/ratsig1
1089 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1090 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1091 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1095 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1096 sigmaii(i,j)=rsum_max
1097 sigmaii(j,i)=rsum_max
1099 c sigmaii(i,j)=r0(i,j)
1100 c sigmaii(j,i)=r0(i,j)
1102 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1103 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1104 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1105 augm(i,j)=epsij*r_augm**(2*expon)
1106 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1113 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1114 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1115 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1121 C Define the SC-p interaction constants (hard-coded; old style)
1124 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1126 c aad(i,1)=0.3D0*4.0D0**12
1127 C Following line for constants currently implemented
1128 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1129 aad(i,1)=1.5D0*4.0D0**12
1130 c aad(i,1)=0.17D0*5.6D0**12
1132 C "Soft" SC-p repulsion
1134 C Following line for constants currently implemented
1135 c aad(i,1)=0.3D0*4.0D0**6
1136 C "Hard" SC-p repulsion
1137 bad(i,1)=3.0D0*4.0D0**6
1138 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1147 C 8/9/01 Read the SC-p interaction constants from file
1150 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1153 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1154 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1155 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1156 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1160 write (iout,*) "Parameters of SC-p interactions:"
1162 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1163 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1168 C Define the constants of the disulfide bridge
1172 c Old arbitrary potential - commented out.
1177 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1178 c energy surface of diethyl disulfide.
1179 c A. Liwo and U. Kozlowska, 11/24/03
1196 write (iout,'(/a)') "Disulfide bridge parameters:"
1197 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1198 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1199 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1200 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1204 111 write (iout,*) "Error reading bending energy parameters."
1206 112 write (iout,*) "Error reading rotamer energy parameters."
1208 113 write (iout,*) "Error reading torsional energy parameters."
1211 & "Error reading side-chain torsional energy parameters."
1213 114 write (iout,*) "Error reading double torsional energy parameters."
1216 & "Error reading cumulant (multibody energy) parameters."
1218 116 write (iout,*) "Error reading electrostatic energy parameters."
1220 117 write (iout,*) "Error reading side chain interaction parameters."
1222 118 write (iout,*) "Error reading SCp interaction parameters."
1224 119 write (iout,*) "Error reading SCCOR parameters"
1227 call MPI_Finalize(Ierror)
1234 subroutine getenv_loc(var, val)
1235 character(*) var, val
1238 character(2000) line
1241 open (196,file='env',status='old',readonly,shared)
1243 c write(*,*)'looking for ',var
1244 10 read(196,*,err=11,end=11)line
1245 iread=index(line,var)
1246 c write(*,*)iread,' ',var,' ',line
1247 if (iread.eq.0) go to 10
1248 c write(*,*)'---> ',line
1254 iread=iread+ilen(var)+1
1255 read (line(iread:),*,err=12,end=12) val
1256 c write(*,*)'OK: ',var,' = ',val
1262 #elif (defined CRAY)
1263 integer lennam,lenval,ierror
1265 c getenv using a POSIX call, useful on the T3D
1266 c Sept 1996, comment out error check on advice of H. Pritchard
1269 if(lennam.le.0) stop '--error calling getenv--'
1270 call pxfgetenv(var,lennam,val,lenval,ierror)
1271 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1273 call getenv(var,val)