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)
151 & 'Parameters of the virtual-bond valence angles:'
152 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
153 & ' ATHETA0 ',' A1 ',' A2 ',
156 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
157 & a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
159 write (iout,'(/a/9x,5a/79(1h-))')
160 & 'Parameters of the expression for sigma(theta_c):',
161 & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
162 & ' ALPH3 ',' SIGMA0C '
164 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
165 & (polthet(j,i),j=0,3),sigc0(i)
167 write (iout,'(/a/9x,5a/79(1h-))')
168 & 'Parameters of the second gaussian:',
169 & ' THETA0 ',' SIGMA0 ',' G1 ',
172 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
173 & sig0(i),(gthet(j,i),j=1,3)
177 & 'Parameters of the virtual-bond valence angles:'
178 write (iout,'(/a/9x,5a/79(1h-))')
179 & 'Coefficients of expansion',
180 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
181 & ' b1*10^1 ',' b2*10^1 '
183 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
184 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
185 $ (10*bthet(j,i,1,1),j=1,2)
187 write (iout,'(/a/9x,5a/79(1h-))')
188 & 'Parameters of the expression for sigma(theta_c):',
189 & ' alpha0 ',' alph1 ',' alph2 ',
190 & ' alhp3 ',' sigma0c '
192 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
193 & (polthet(j,i),j=0,3),sigc0(i)
195 write (iout,'(/a/9x,5a/79(1h-))')
196 & 'Parameters of the second gaussian:',
197 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
200 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
201 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
207 C Read the parameters of Utheta determined from ab initio surfaces
208 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
210 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,
211 & ntheterm3,nsingle,ndouble
212 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
213 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
219 aathet(l,i,j,k)=0.0d0
223 bbthet(m,l,i,j,k)=0.0d0
224 ccthet(m,l,i,j,k)=0.0d0
225 ddthet(m,l,i,j,k)=0.0d0
226 eethet(m,l,i,j,k)=0.0d0
232 ffthet(mm,m,l,i,j,k)=0.0d0
233 ggthet(mm,m,l,i,j,k)=0.0d0
243 read (ithep,'(3a)',end=111,err=111) res1,res2,res3
244 read (ithep,*,end=111,err=111) aa0thet(i,j,k)
245 read (ithep,*,end=111,err=111)(aathet(l,i,j,k),l=1,ntheterm)
246 read (ithep,*,end=111,err=111)
247 & ((bbthet(lll,ll,i,j,k),lll=1,nsingle),
248 & (ccthet(lll,ll,i,j,k),lll=1,nsingle),
249 & (ddthet(lll,ll,i,j,k),lll=1,nsingle),
250 & (eethet(lll,ll,i,j,k),lll=1,nsingle),ll=1,ntheterm2)
251 read (ithep,*,end=111,err=111)
252 & (((ffthet(llll,lll,ll,i,j,k),ffthet(lll,llll,ll,i,j,k),
253 & ggthet(llll,lll,ll,i,j,k),ggthet(lll,llll,ll,i,j,k),
254 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
259 C For dummy ends assign glycine-type coefficients of theta-only terms; the
260 C coefficients of theta-and-gamma-dependent terms are zero.
265 aathet(l,i,j,nthetyp+1)=aathet(l,i,j,1)
266 aathet(l,nthetyp+1,i,j)=aathet(l,1,i,j)
268 aa0thet(i,j,nthetyp+1)=aa0thet(i,j,1)
269 aa0thet(nthetyp+1,i,j)=aa0thet(1,i,j)
272 aathet(l,nthetyp+1,i,nthetyp+1)=aathet(l,1,i,1)
274 aa0thet(nthetyp+1,i,nthetyp+1)=aa0thet(1,i,1)
277 C Control printout of the coefficients of virtual-bond-angle potentials
280 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
284 write (iout,'(//4a)')
285 & 'Type ',onelett(i),onelett(j),onelett(k)
286 write (iout,'(//a,10x,a)') " l","a[l]"
287 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k)
288 write (iout,'(i2,1pe15.5)')
289 & (l,aathet(l,i,j,k),l=1,ntheterm)
291 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))')
292 & "b",l,"c",l,"d",l,"e",l
294 write (iout,'(i2,4(1pe15.5))') m,
295 & bbthet(m,l,i,j,k),ccthet(m,l,i,j,k),
296 & ddthet(m,l,i,j,k),eethet(m,l,i,j,k)
300 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
301 & "f+",l,"f-",l,"g+",l,"g-",l
304 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
305 & ffthet(n,m,l,i,j,k),ffthet(m,n,l,i,j,k),
306 & ggthet(n,m,l,i,j,k),ggthet(m,n,l,i,j,k)
315 write (2,*) "Start reading THETA_PDB"
317 read (ithep_pdb,*,err=111,end=111) a0thet(i),(athet(j,i),j=1,2),
319 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
320 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
321 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
323 write (2,*) "End reading THETA_PDB"
330 C Read the parameters of the probability distribution/energy expression
331 C of the side chains.
334 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
338 dsc_inv(i)=1.0D0/dsc(i)
349 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
350 & ((blower(k,l,1),l=1,k),k=1,3)
351 censc(1,1,-i)=censc(1,1,i)
352 censc(2,1,-i)=censc(2,1,i)
353 censc(3,1,-i)=-censc(3,1,i)
356 read (irotam,*,end=112,err=112) bsc(j,i)
357 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
358 & ((blower(k,l,j),l=1,k),k=1,3)
359 censc(1,j,-i)=censc(1,j,i)
360 censc(2,j,-i)=censc(2,j,i)
361 censc(3,j,-i)=-censc(3,j,i)
362 C BSC is amplitude of Gaussian
369 akl=akl+blower(k,m,j)*blower(l,m,j)
373 if (((k.eq.3).and.(l.ne.3))
374 & .or.((l.eq.3).and.(k.ne.3))) then
375 gaussc(k,l,j,-i)=-akl
376 gaussc(l,k,j,-i)=-akl
388 write (iout,'(/a)') 'Parameters of side-chain local geometry'
393 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
394 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
395 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
396 & 'log h',(bsc(j,i),j=1,nlobi)
397 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
398 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
400 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
401 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
404 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
405 write (iout,'(a,f10.4,4(16x,f10.4))')
406 & 'Center ',(bsc(j,i),j=1,nlobi)
407 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
416 C Read scrot parameters for potentials determined from all-atom AM1 calculations
417 C added by Urszula Kozlowska 07/11/2007
420 read (irotam,*,end=112,err=112)
422 read (irotam,*,end=112,err=112)
425 read(irotam,*,end=112,err=112) sc_parmin(j,i)
430 C Read the parameters of the probability distribution/energy expression
431 C of the side chains.
434 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
438 dsc_inv(i)=1.0D0/dsc(i)
449 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
450 & ((blower(k,l,1),l=1,k),k=1,3)
452 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
453 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
454 & ((blower(k,l,j),l=1,k),k=1,3)
461 akl=akl+blower(k,m,j)*blower(l,m,j)
476 C Read torsional parameters in old format
478 read (itorp,*,end=113,err=113) ntortyp,nterm_old
479 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
480 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
485 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
491 write (iout,'(/a/)') 'Torsional constants:'
494 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
495 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
501 C Read torsional parameters
503 read (itorp,*,end=113,err=113) ntortyp
504 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
508 itortyp(i)=-itortyp(-1)
510 c write (iout,*) 'ntortyp',ntortyp
511 do i=0,ntortyp,ntortyp-1
512 do j=-ntortyp,ntortyp
513 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
515 nterm(-i,-j,iblock)=nterm(i,j,iblock)
516 nlor(-i,-j,iblock)=nlor(i,j,iblock)
518 itortyp(i)=-itortyp(-i)
520 c write (iout,*) 'ntortyp',ntortyp
522 do j=-ntortyp,ntortyp
523 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
525 nterm(-i,-j,iblock)=nterm(i,j,iblock)
530 do k=1,nterm(i,j,iblock)
531 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
533 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
535 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
537 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
539 v0ij=v0ij+si*v1(k,i,j,iblock)
541 c write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
542 c write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
543 c &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
545 do k=1,nlor(i,j,iblock)
546 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
547 & vlor2(k,i,j),vlor3(k,i,j)
548 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
551 v0(-i,-j,iblock)=v0ij
557 write (iout,'(/a/)') 'Torsional constants:'
560 write (iout,*) 'ityp',i,' jtyp',j
561 write (iout,*) 'Fourier constants'
562 do k=1,nterm(i,j,iblock)
563 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
566 write (iout,*) 'Lorenz constants'
567 do k=1,nlor(i,j,iblock)
568 write (iout,'(3(1pe15.5))')
569 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
575 C 6/23/01 Read parameters for double torsionals
579 do j=-ntortyp+1,ntortyp-1
580 do k=-ntortyp+1,ntortyp-1
581 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
582 write (iout,*) "OK onelett",
585 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
586 & .or. t3.ne.toronelet(k)) then
587 write (iout,*) "Error in double torsional parameter file",
590 call MPI_Finalize(Ierror)
592 stop "Error in double torsional parameter file"
594 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
595 & ntermd_2(i,j,k,iblock)
596 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
597 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
598 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
599 & ntermd_1(i,j,k,iblock))
600 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
601 & ntermd_1(i,j,k,iblock))
602 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
603 & ntermd_1(i,j,k,iblock))
604 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
605 & ntermd_1(i,j,k,iblock))
606 C Matrix of D parameters for one dimesional foureir series
607 do l=1,ntermd_1(i,j,k,iblock)
608 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
609 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
610 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
611 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
612 c write(iout,*) "whcodze" ,
613 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
615 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
616 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
617 & v2s(m,l,i,j,k,iblock),
618 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
619 C Matrix of D parameters for two dimesional fourier series
620 do l=1,ntermd_2(i,j,k,iblock)
622 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
623 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
624 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
625 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
634 write (iout,*) 'Constants for double torsionals'
637 do j=-ntortyp,ntortyp
638 do k=-ntortyp,ntortyp
639 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
640 & ' nsingle',ntermd_1(i,j,k,iblock),
641 & ' ndouble',ntermd_2(i,j,k,iblock)
643 write (iout,*) 'Single angles:'
644 do l=1,ntermd_1(i,j,k,iblock)
645 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
646 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
647 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock)
650 write (iout,*) 'Pairs of angles:'
651 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
652 do l=1,ntermd_2(i,j,k,iblock)
653 write (iout,'(i5,20f10.5)')
654 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
657 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
658 do l=1,ntermd_2(i,j,k,iblock)
659 write (iout,'(i5,20f10.5)')
660 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
661 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
671 C 5/21/07 (AL) Read coefficients of the backbone-local sidechain-local
672 C correlation energies.
674 read (isccor,*,end=119,err=119) nterm_sccor
679 read (isccor,*,end=119,err=119) kk,v1sccor(k,i,j),
686 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
689 write (iout,*) 'ityp',i,' jtyp',j
691 write (iout,'(2(1pe15.5))') v1sccor(k,i,j),v2sccor(k,i,j)
697 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
698 C interaction energy of the Gly, Ala, and Pro prototypes.
702 write (iout,*) "Coefficients of the cumulants"
704 read (ifourier,*) nloctyp
707 read (ifourier,*,end=115,err=115)
708 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
710 write (iout,*) 'Type',i
711 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
752 c Ctilde(1,1,i)=0.0d0
753 c Ctilde(1,2,i)=0.0d0
754 c Ctilde(2,1,i)=0.0d0
755 c Ctilde(2,2,i)=0.0d0
776 c Dtilde(1,1,i)=0.0d0
777 c Dtilde(1,2,i)=0.0d0
778 c Dtilde(2,1,i)=0.0d0
779 c Dtilde(2,2,i)=0.0d0
780 EE(1,1,i)= b(10)+b(11)
781 EE(2,2,i)=-b(10)+b(11)
782 EE(2,1,i)= b(12)-b(13)
783 EE(1,2,i)= b(12)+b(13)
784 EE(1,1,-i)= b(10)+b(11)
785 EE(2,2,-i)=-b(10)+b(11)
786 EE(2,1,-i)=-b(12)+b(13)
787 EE(1,2,-i)=-b(12)-b(13)
792 c ee(2,1,i)=ee(1,2,i)
796 write (iout,*) 'Type',i
798 write(iout,*) B1(1,i),B1(2,i)
800 write(iout,*) B2(1,i),B2(2,i)
803 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
807 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
811 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
816 C Read electrostatic-interaction parameters
820 write (iout,'(/a)') 'Electrostatic interaction constants:'
821 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
822 & 'IT','JT','APP','BPP','AEL6','AEL3'
824 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
825 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
826 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
827 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
832 app (i,j)=epp(i,j)*rri*rri
833 bpp (i,j)=-2.0D0*epp(i,j)*rri
834 ael6(i,j)=elpp6(i,j)*4.2D0**6
835 ael3(i,j)=elpp3(i,j)*4.2D0**3
836 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
837 & ael6(i,j),ael3(i,j)
841 C Read side-chain interaction parameters.
843 read (isidep,*,end=117,err=117) ipot,expon
844 if (ipot.lt.1 .or. ipot.gt.5) then
845 write (iout,'(2a)') 'Error while reading SC interaction',
846 & 'potential file - unknown potential type.'
848 call MPI_Finalize(Ierror)
854 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
855 & ', exponents are ',expon,2*expon
856 goto (10,20,30,30,40) ipot
857 C----------------------- LJ potential ---------------------------------
858 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
859 & (sigma0(i),i=1,ntyp)
861 write (iout,'(/a/)') 'Parameters of the LJ potential:'
862 write (iout,'(a/)') 'The epsilon array:'
863 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
864 write (iout,'(/a)') 'One-body parameters:'
865 write (iout,'(a,4x,a)') 'residue','sigma'
866 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
869 C----------------------- LJK potential --------------------------------
870 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
871 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
873 write (iout,'(/a/)') 'Parameters of the LJK potential:'
874 write (iout,'(a/)') 'The epsilon array:'
875 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
876 write (iout,'(/a)') 'One-body parameters:'
877 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
878 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
882 C---------------------- GB or BP potential -----------------------------
883 30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
884 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
886 C For the GB potential convert sigma'**2 into chi'
889 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
893 write (iout,'(/a/)') 'Parameters of the BP potential:'
894 write (iout,'(a/)') 'The epsilon array:'
895 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
896 write (iout,'(/a)') 'One-body parameters:'
897 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
899 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
900 & chip(i),alp(i),i=1,ntyp)
903 C--------------------- GBV potential -----------------------------------
904 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
905 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
906 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
908 write (iout,'(/a/)') 'Parameters of the GBV potential:'
909 write (iout,'(a/)') 'The epsilon array:'
910 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
911 write (iout,'(/a)') 'One-body parameters:'
912 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
913 & 's||/s_|_^2',' chip ',' alph '
914 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
915 & sigii(i),chip(i),alp(i),i=1,ntyp)
919 C-----------------------------------------------------------------------
920 C Calculate the "working" parameters of SC interactions.
928 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
929 sigma(j,i)=sigma(i,j)
930 rs0(i,j)=dwa16*sigma(i,j)
934 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
935 & 'Working parameters of the SC interactions:',
936 & ' a ',' b ',' augm ',' sigma ',' r0 ',
941 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
950 sigeps=dsign(1.0D0,epsij)
952 aa(i,j)=epsij*rrij*rrij
953 bb(i,j)=-sigeps*epsij*rrij
961 ratsig1=sigt2sq/sigt1sq
962 ratsig2=1.0D0/ratsig1
963 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
964 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
965 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
969 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
970 sigmaii(i,j)=rsum_max
971 sigmaii(j,i)=rsum_max
973 c sigmaii(i,j)=r0(i,j)
974 c sigmaii(j,i)=r0(i,j)
976 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
977 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
978 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
979 augm(i,j)=epsij*r_augm**(2*expon)
980 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
987 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
988 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
989 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
995 C Define the SC-p interaction constants (hard-coded; old style)
998 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1000 c aad(i,1)=0.3D0*4.0D0**12
1001 C Following line for constants currently implemented
1002 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1003 aad(i,1)=1.5D0*4.0D0**12
1004 c aad(i,1)=0.17D0*5.6D0**12
1006 C "Soft" SC-p repulsion
1008 C Following line for constants currently implemented
1009 c aad(i,1)=0.3D0*4.0D0**6
1010 C "Hard" SC-p repulsion
1011 bad(i,1)=3.0D0*4.0D0**6
1012 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1021 C 8/9/01 Read the SC-p interaction constants from file
1024 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1027 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1028 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1029 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1030 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1034 write (iout,*) "Parameters of SC-p interactions:"
1036 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1037 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1042 C Define the constants of the disulfide bridge
1046 c Old arbitrary potential - commented out.
1051 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1052 c energy surface of diethyl disulfide.
1053 c A. Liwo and U. Kozlowska, 11/24/03
1070 write (iout,'(/a)') "Disulfide bridge parameters:"
1071 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1072 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1073 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1074 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1078 111 write (iout,*) "Error reading bending energy parameters."
1080 112 write (iout,*) "Error reading rotamer energy parameters."
1082 113 write (iout,*) "Error reading torsional energy parameters."
1084 114 write (iout,*) "Error reading double torsional energy parameters."
1087 & "Error reading cumulant (multibody energy) parameters."
1089 116 write (iout,*) "Error reading electrostatic energy parameters."
1091 117 write (iout,*) "Error reading side chain interaction parameters."
1093 118 write (iout,*) "Error reading SCp interaction parameters."
1095 119 write (iout,*) "Error reading SCCOR parameters"
1098 call MPI_Finalize(Ierror)
1105 subroutine getenv_loc(var, val)
1106 character(*) var, val
1109 character(2000) line
1112 open (196,file='env',status='old',readonly,shared)
1114 c write(*,*)'looking for ',var
1115 10 read(196,*,err=11,end=11)line
1116 iread=index(line,var)
1117 c write(*,*)iread,' ',var,' ',line
1118 if (iread.eq.0) go to 10
1119 c write(*,*)'---> ',line
1125 iread=iread+ilen(var)+1
1126 read (line(iread:),*,err=12,end=12) val
1127 c write(*,*)'OK: ',var,' = ',val
1133 #elif (defined CRAY)
1134 integer lennam,lenval,ierror
1136 c getenv using a POSIX call, useful on the T3D
1137 c Sept 1996, comment out error check on advice of H. Pritchard
1140 if(lennam.le.0) stop '--error calling getenv--'
1141 call pxfgetenv(var,lennam,val,lenval,ierror)
1142 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1144 call getenv(var,val)