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)
324 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)
507 itortyp(i)=-itortyp(-i)
509 c write (iout,*) 'ntortyp',ntortyp
511 do j=-ntortyp+1,ntortyp-1
512 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
514 nterm(-i,-j,iblock)=nterm(i,j,iblock)
515 nlor(-i,-j,iblock)=nlor(i,j,iblock)
518 do k=1,nterm(i,j,iblock)
519 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
521 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
522 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
523 v0ij=v0ij+si*v1(k,i,j,iblock)
525 c write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
526 c write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
527 c &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
529 do k=1,nlor(i,j,iblock)
530 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
531 & vlor2(k,i,j),vlor3(k,i,j)
532 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
535 v0(-i,-j,iblock)=v0ij
541 write (iout,'(/a/)') 'Torsional constants:'
544 write (iout,*) 'ityp',i,' jtyp',j
545 write (iout,*) 'Fourier constants'
546 do k=1,nterm(i,j,iblock)
547 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
550 write (iout,*) 'Lorenz constants'
551 do k=1,nlor(i,j,iblock)
552 write (iout,'(3(1pe15.5))')
553 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
559 C 6/23/01 Read parameters for double torsionals
563 do j=-ntortyp+1,ntortyp-1
564 do k=-ntortyp+1,ntortyp-1
565 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
566 c write (iout,*) "OK onelett",
569 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
570 & .or. t3.ne.toronelet(k)) then
571 write (iout,*) "Error in double torsional parameter file",
574 call MPI_Finalize(Ierror)
576 stop "Error in double torsional parameter file"
578 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
579 & ntermd_2(i,j,k,iblock)
580 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
581 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
582 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
583 & ntermd_1(i,j,k,iblock))
584 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
585 & ntermd_1(i,j,k,iblock))
586 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
587 & ntermd_1(i,j,k,iblock))
588 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
589 & ntermd_1(i,j,k,iblock))
590 C Martix of D parameters for one dimesional foureir series
591 do l=1,ntermd_1(i,j,k,iblock)
592 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
593 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
594 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
595 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
596 c write(iout,*) "whcodze" ,
597 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
599 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
600 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
601 & v2s(m,l,i,j,k,iblock),
602 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
603 C Martix of D parameters for two dimesional fourier series
604 do l=1,ntermd_2(i,j,k,iblock)
606 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
607 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
608 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
609 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
618 write (iout,*) 'Constants for double torsionals'
621 do j=-ntortyp,ntortyp
622 do k=-ntortyp,ntortyp
623 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
624 & ' nsingle',ntermd_1(i,j,k,iblock),
625 & ' ndouble',ntermd_2(i,j,k,iblock)
627 write (iout,*) 'Single angles:'
628 do l=1,ntermd_1(i,j,k,iblock)
629 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
630 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
631 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
632 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
635 write (iout,*) 'Pairs of angles:'
636 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
637 do l=1,ntermd_2(i,j,k,iblock)
638 write (iout,'(i5,20f10.5)')
639 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
642 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
643 do l=1,ntermd_2(i,j,k,iblock)
644 write (iout,'(i5,20f10.5)')
645 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
646 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
655 C Read of Side-chain backbone correlation parameters
656 C Modified 11 May 2012 by Adasko
659 read (isccor,*,end=113,err=113) nsccortyp
661 read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
663 isccortyp(i)=-isccortyp(-i)
665 iscprol=isccortyp(20)
666 c write (iout,*) 'ntortyp',ntortyp
668 cc maxinter is maximum interaction sites
672 read (isccor,*,end=113,err=113) nterm_sccor(i,j),nlor_sccor(i,j)
678 nterm_sccor(-i,j)=nterm_sccor(i,j)
679 nterm_sccor(-i,-j)=nterm_sccor(i,j)
680 nterm_sccor(i,-j)=nterm_sccor(i,j)
681 do k=1,nterm_sccor(i,j)
682 read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j)
684 if (j.eq.iscprol) then
685 if (i.eq.isccortyp(10)) then
686 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
687 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
689 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
690 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
691 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
692 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
693 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
694 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
695 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
696 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
699 if (i.eq.isccortyp(10)) then
700 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
701 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
703 if (j.eq.isccortyp(10)) then
704 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
705 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
707 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
708 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
709 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
710 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
711 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
712 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
716 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
717 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
718 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
719 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
722 do k=1,nlor_sccor(i,j)
723 read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
724 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
725 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
726 &(1+vlor3sccor(k,i,j)**2)
728 v0sccor(l,i,j)=v0ijsccor
729 v0sccor(l,-i,j)=v0ijsccor1
730 v0sccor(l,i,-j)=v0ijsccor2
731 v0sccor(l,-i,-j)=v0ijsccor3
737 read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
738 c write (iout,*) 'ntortyp',ntortyp
740 cc maxinter is maximum interaction sites
744 read (isccor,*,end=113,err=113)
745 & nterm_sccor(i,j),nlor_sccor(i,j)
749 do k=1,nterm_sccor(i,j)
750 read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j)
752 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
755 do k=1,nlor_sccor(i,j)
756 read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
757 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
758 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
759 &(1+vlor3sccor(k,i,j)**2)
761 v0sccor(i,j)=v0ijsccor
769 write (iout,'(/a/)') 'Torsional constants:'
772 write (iout,*) 'ityp',i,' jtyp',j
773 write (iout,*) 'Fourier constants'
774 do k=1,nterm_sccor(i,j)
775 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
777 write (iout,*) 'Lorenz constants'
778 do k=1,nlor_sccor(i,j)
779 write (iout,'(3(1pe15.5))')
780 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
787 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
788 C interaction energy of the Gly, Ala, and Pro prototypes.
792 write (iout,*) "Coefficients of the cumulants"
794 read (ifourier,*) nloctyp
796 read (ifourier,*,end=115,err=115)
797 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
799 write (iout,*) 'Type',i
800 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
842 c Ctilde(1,1,i)=0.0d0
843 c Ctilde(1,2,i)=0.0d0
844 c Ctilde(2,1,i)=0.0d0
845 c Ctilde(2,2,i)=0.0d0
867 c Dtilde(1,1,i)=0.0d0
868 c Dtilde(1,2,i)=0.0d0
869 c Dtilde(2,1,i)=0.0d0
870 c Dtilde(2,2,i)=0.0d0
871 EE(1,1,i)= b(10)+b(11)
872 EE(2,2,i)=-b(10)+b(11)
873 EE(2,1,i)= b(12)-b(13)
874 EE(1,2,i)= b(12)+b(13)
875 EE(1,1,-i)= b(10)+b(11)
876 EE(2,2,-i)=-b(10)+b(11)
877 EE(2,1,-i)=-b(12)+b(13)
878 EE(1,2,-i)=-b(12)-b(13)
884 c ee(2,1,i)=ee(1,2,i)
888 write (iout,*) 'Type',i
890 write(iout,*) B1(1,i),B1(2,i)
892 write(iout,*) B2(1,i),B2(2,i)
895 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
899 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
903 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
908 C Read electrostatic-interaction parameters
912 write (iout,'(/a)') 'Electrostatic interaction constants:'
913 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
914 & 'IT','JT','APP','BPP','AEL6','AEL3'
916 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
917 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
918 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
919 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
924 app (i,j)=epp(i,j)*rri*rri
925 bpp (i,j)=-2.0D0*epp(i,j)*rri
926 ael6(i,j)=elpp6(i,j)*4.2D0**6
927 ael3(i,j)=elpp3(i,j)*4.2D0**3
928 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
929 & ael6(i,j),ael3(i,j)
933 C Read side-chain interaction parameters.
935 read (isidep,*,end=117,err=117) ipot,expon
936 if (ipot.lt.1 .or. ipot.gt.5) then
937 write (iout,'(2a)') 'Error while reading SC interaction',
938 & 'potential file - unknown potential type.'
940 call MPI_Finalize(Ierror)
946 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
947 & ', exponents are ',expon,2*expon
948 goto (10,20,30,30,40) ipot
949 C----------------------- LJ potential ---------------------------------
950 10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
951 & (sigma0(i),i=1,ntyp)
953 write (iout,'(/a/)') 'Parameters of the LJ potential:'
954 write (iout,'(a/)') 'The epsilon array:'
955 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
956 write (iout,'(/a)') 'One-body parameters:'
957 write (iout,'(a,4x,a)') 'residue','sigma'
958 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
961 C----------------------- LJK potential --------------------------------
962 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
963 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
965 write (iout,'(/a/)') 'Parameters of the LJK potential:'
966 write (iout,'(a/)') 'The epsilon array:'
967 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
968 write (iout,'(/a)') 'One-body parameters:'
969 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
970 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
974 C---------------------- GB or BP potential -----------------------------
975 30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
976 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
978 C For the GB potential convert sigma'**2 into chi'
981 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
985 write (iout,'(/a/)') 'Parameters of the BP potential:'
986 write (iout,'(a/)') 'The epsilon array:'
987 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
988 write (iout,'(/a)') 'One-body parameters:'
989 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
991 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
992 & chip(i),alp(i),i=1,ntyp)
995 C--------------------- GBV potential -----------------------------------
996 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
997 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
998 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1000 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1001 write (iout,'(a/)') 'The epsilon array:'
1002 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1003 write (iout,'(/a)') 'One-body parameters:'
1004 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1005 & 's||/s_|_^2',' chip ',' alph '
1006 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1007 & sigii(i),chip(i),alp(i),i=1,ntyp)
1011 C-----------------------------------------------------------------------
1012 C Calculate the "working" parameters of SC interactions.
1020 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1021 sigma(j,i)=sigma(i,j)
1022 rs0(i,j)=dwa16*sigma(i,j)
1026 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1027 & 'Working parameters of the SC interactions:',
1028 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1033 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1042 sigeps=dsign(1.0D0,epsij)
1044 aa(i,j)=epsij*rrij*rrij
1045 bb(i,j)=-sigeps*epsij*rrij
1049 sigt1sq=sigma0(i)**2
1050 sigt2sq=sigma0(j)**2
1053 ratsig1=sigt2sq/sigt1sq
1054 ratsig2=1.0D0/ratsig1
1055 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1056 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1057 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1061 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1062 sigmaii(i,j)=rsum_max
1063 sigmaii(j,i)=rsum_max
1065 c sigmaii(i,j)=r0(i,j)
1066 c sigmaii(j,i)=r0(i,j)
1068 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1069 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1070 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1071 augm(i,j)=epsij*r_augm**(2*expon)
1072 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1079 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1080 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1081 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1087 C Define the SC-p interaction constants (hard-coded; old style)
1090 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1092 c aad(i,1)=0.3D0*4.0D0**12
1093 C Following line for constants currently implemented
1094 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1095 aad(i,1)=1.5D0*4.0D0**12
1096 c aad(i,1)=0.17D0*5.6D0**12
1098 C "Soft" SC-p repulsion
1100 C Following line for constants currently implemented
1101 c aad(i,1)=0.3D0*4.0D0**6
1102 C "Hard" SC-p repulsion
1103 bad(i,1)=3.0D0*4.0D0**6
1104 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1113 C 8/9/01 Read the SC-p interaction constants from file
1116 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1119 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1120 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1121 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1122 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1126 write (iout,*) "Parameters of SC-p interactions:"
1128 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1129 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1134 C Define the constants of the disulfide bridge
1138 c Old arbitrary potential - commented out.
1143 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1144 c energy surface of diethyl disulfide.
1145 c A. Liwo and U. Kozlowska, 11/24/03
1162 write (iout,'(/a)') "Disulfide bridge parameters:"
1163 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1164 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1165 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1166 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1170 111 write (iout,*) "Error reading bending energy parameters."
1172 112 write (iout,*) "Error reading rotamer energy parameters."
1174 113 write (iout,*) "Error reading torsional energy parameters."
1176 114 write (iout,*) "Error reading double torsional energy parameters."
1179 & "Error reading cumulant (multibody energy) parameters."
1181 116 write (iout,*) "Error reading electrostatic energy parameters."
1183 117 write (iout,*) "Error reading side chain interaction parameters."
1185 118 write (iout,*) "Error reading SCp interaction parameters."
1187 119 write (iout,*) "Error reading SCCOR parameters"
1190 call MPI_Finalize(Ierror)
1197 subroutine getenv_loc(var, val)
1198 character(*) var, val
1201 character(2000) line
1204 open (196,file='env',status='old',readonly,shared)
1206 c write(*,*)'looking for ',var
1207 10 read(196,*,err=11,end=11)line
1208 iread=index(line,var)
1209 c write(*,*)iread,' ',var,' ',line
1210 if (iread.eq.0) go to 10
1211 c write(*,*)'---> ',line
1217 iread=iread+ilen(var)+1
1218 read (line(iread:),*,err=12,end=12) val
1219 c write(*,*)'OK: ',var,' = ',val
1225 #elif (defined CRAY)
1226 integer lennam,lenval,ierror
1228 c getenv using a POSIX call, useful on the T3D
1229 c Sept 1996, comment out error check on advice of H. Pritchard
1232 if(lennam.le.0) stop '--error calling getenv--'
1233 call pxfgetenv(var,lennam,val,lenval,ierror)
1234 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1236 call getenv(var,val)