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 implicit real*8 (a-h,o-z)
9 include 'COMMON.IOUNITS'
10 include 'COMMON.CHAIN'
11 include 'COMMON.INTERACT'
13 include 'COMMON.LOCAL'
14 include 'COMMON.TORSION'
15 include 'COMMON.FFIELD'
16 include 'COMMON.NAMES'
17 include 'COMMON.SBRIDGE'
18 include 'COMMON.SCCOR'
19 include 'COMMON.SCROT'
21 character*1 onelett(-2:2) /"p","a","G","A","P"/
23 dimension blower(3,3,maxlob)
24 double precision ip,mp
28 C Set LPRINT=.TRUE. for debugging
29 dwa16=2.0d0**(1.0d0/6.0d0)
32 C Assign virtual-bond length
37 read (ibond,*) vbldp0,akp
40 read (ibond,*) vbldsc0(1,i),aksc(1,i)
45 dsc_inv(i)=1.0D0/dsc(i)
49 read (ibond,*) ijunk,vbldp0,akp,rjunk
51 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
57 dsc_inv(i)=1.0D0/dsc(i)
62 write(iout,'(/a/)')"Force constants virtual bonds:"
63 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
65 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
67 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
68 & vbldsc0(1,i),aksc(1,i),abond0(1,i)
70 write (iout,'(13x,3f10.5)')
71 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
77 C Read the parameters of the probability distribution/energy expression
78 C of the virtual-bond valence angles theta
81 read (ithep,*) a0thet(i),(athet(j,i,1,1),j=1,2),
82 & (bthet(j,i,1,1),j=1,2)
83 read (ithep,*) (polthet(j,i),j=0,3)
84 read (ithep,*) (gthet(j,i),j=1,3)
85 read (ithep,*) theta0(i),sig0(i),sigc0(i)
89 athet(1,i,1,-1)=athet(1,i,1,1)
90 athet(2,i,1,-1)=athet(2,i,1,1)
91 bthet(1,i,1,-1)=-bthet(1,i,1,1)
92 bthet(2,i,1,-1)=-bthet(2,i,1,1)
93 athet(1,i,-1,1)=-athet(1,i,1,1)
94 athet(2,i,-1,1)=-athet(2,i,1,1)
95 bthet(1,i,-1,1)=bthet(1,i,1,1)
96 bthet(2,i,-1,1)=bthet(2,i,1,1)
100 athet(1,i,-1,-1)=athet(1,-i,1,1)
101 athet(2,i,-1,-1)=-athet(2,-i,1,1)
102 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
103 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
104 athet(1,i,-1,1)=athet(1,-i,1,1)
105 athet(2,i,-1,1)=-athet(2,-i,1,1)
106 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
107 bthet(2,i,-1,1)=bthet(2,-i,1,1)
108 athet(1,i,1,-1)=-athet(1,-i,1,1)
109 athet(2,i,1,-1)=athet(2,-i,1,1)
110 bthet(1,i,1,-1)=bthet(1,-i,1,1)
111 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
116 polthet(j,i)=polthet(j,-i)
119 gthet(j,i)=gthet(j,-i)
125 c & 'Parameters of the virtual-bond valence angles:'
126 c write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
127 c & ' ATHETA0 ',' A1 ',' A2 ',
130 c write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
131 c & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
133 c write (iout,'(/a/9x,5a/79(1h-))')
134 c & 'Parameters of the expression for sigma(theta_c):',
135 c & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
136 c & ' ALPH3 ',' SIGMA0C '
138 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
139 c & (polthet(j,i),j=0,3),sigc0(i)
141 c write (iout,'(/a/9x,5a/79(1h-))')
142 c & 'Parameters of the second gaussian:',
143 c & ' THETA0 ',' SIGMA0 ',' G1 ',
146 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
147 c & sig0(i),(gthet(j,i),j=1,3)
150 & 'Parameters of the virtual-bond valence angles:'
151 write (iout,'(/a/9x,5a/79(1h-))')
152 & 'Coefficients of expansion',
153 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
154 & ' b1*10^1 ',' b2*10^1 '
156 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
157 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
158 & (10*bthet(j,i,1,1),j=1,2)
160 write (iout,'(/a/9x,5a/79(1h-))')
161 & 'Parameters of the expression for sigma(theta_c):',
162 & ' alpha0 ',' alph1 ',' alph2 ',
163 & ' alhp3 ',' sigma0c '
165 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
166 & (polthet(j,i),j=0,3),sigc0(i)
168 write (iout,'(/a/9x,5a/79(1h-))')
169 & 'Parameters of the second gaussian:',
170 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
173 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
174 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
179 C Read the parameters of Utheta determined from ab initio surfaces
180 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
182 read (ithep,*) nthetyp,ntheterm,ntheterm2,
183 & ntheterm3,nsingle,ndouble
184 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
185 read (ithep,*) (ithetyp(i),i=1,ntyp1)
191 aathet(l,i,j,k)=0.0d0
195 bbthet(m,l,i,j,k)=0.0d0
196 ccthet(m,l,i,j,k)=0.0d0
197 ddthet(m,l,i,j,k)=0.0d0
198 eethet(m,l,i,j,k)=0.0d0
204 ffthet(mm,m,l,i,j,k)=0.0d0
205 ggthet(mm,m,l,i,j,k)=0.0d0
215 read (ithep,'(3a)') res1,res2,res3
216 read (ithep,*) aa0thet(i,j,k)
217 read (ithep,*)(aathet(l,i,j,k),l=1,ntheterm)
219 & ((bbthet(lll,ll,i,j,k),lll=1,nsingle),
220 & (ccthet(lll,ll,i,j,k),lll=1,nsingle),
221 & (ddthet(lll,ll,i,j,k),lll=1,nsingle),
222 & (eethet(lll,ll,i,j,k),lll=1,nsingle),ll=1,ntheterm2)
224 & (((ffthet(llll,lll,ll,i,j,k),ffthet(lll,llll,ll,i,j,k),
225 & ggthet(llll,lll,ll,i,j,k),ggthet(lll,llll,ll,i,j,k),
226 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
231 C For dummy ends assign glycine-type coefficients of theta-only terms; the
232 C coefficients of theta-and-gamma-dependent terms are zero.
237 aathet(l,i,j,nthetyp+1)=aathet(l,i,j,1)
238 aathet(l,nthetyp+1,i,j)=aathet(l,1,i,j)
240 aa0thet(i,j,nthetyp+1)=aa0thet(i,j,1)
241 aa0thet(nthetyp+1,i,j)=aa0thet(1,i,j)
244 aathet(l,nthetyp+1,i,nthetyp+1)=aathet(l,1,i,1)
246 aa0thet(nthetyp+1,i,nthetyp+1)=aa0thet(1,i,1)
249 C Control printout of the coefficients of virtual-bond-angle potentials
252 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
256 write (iout,'(//4a)')
257 & 'Type ',onelett(i),onelett(j),onelett(k)
258 write (iout,'(//a,10x,a)') " l","a[l]"
259 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k)
260 write (iout,'(i2,1pe15.5)')
261 & (l,aathet(l,i,j,k),l=1,ntheterm)
263 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
264 & "b",l,"c",l,"d",l,"e",l
266 write (iout,'(i2,4(1pe15.5))') m,
267 & bbthet(m,l,i,j,k),ccthet(m,l,i,j,k),
268 & ddthet(m,l,i,j,k),eethet(m,l,i,j,k)
272 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
273 & "f+",l,"f-",l,"g+",l,"g-",l
276 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
277 & ffthet(n,m,l,i,j,k),ffthet(m,n,l,i,j,k),
278 & ggthet(n,m,l,i,j,k),ggthet(m,n,l,i,j,k)
291 C Read the parameters of the probability distribution/energy expression
292 C of the side chains.
295 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
299 dsc_inv(i)=1.0D0/dsc(i)
310 read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
311 censc(1,1,-i)=censc(1,1,i)
312 censc(2,1,-i)=censc(2,1,i)
313 censc(3,1,-i)=-censc(3,1,i)
315 read (irotam,*) bsc(j,i)
316 read (irotam,*) (censc(k,j,i),k=1,3),
317 & ((blower(k,l,j),l=1,k),k=1,3)
318 censc(1,j,-i)=censc(1,j,i)
319 censc(2,j,-i)=censc(2,j,i)
320 censc(3,j,-i)=-censc(3,j,i)
321 C BSC is amplitude of Gaussian
328 akl=akl+blower(k,m,j)*blower(l,m,j)
332 if (((k.eq.3).and.(l.ne.3))
333 & .or.((l.eq.3).and.(k.ne.3))) then
334 gaussc(k,l,j,-i)=-akl
335 gaussc(l,k,j,-i)=-akl
347 write (iout,'(/a)') 'Parameters of side-chain local geometry'
351 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
352 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
353 c write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
354 c write (iout,'(a,f10.4,4(16x,f10.4))')
355 c & 'Center ',(bsc(j,i),j=1,nlobi)
356 c write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
357 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
358 & 'log h',(bsc(j,i),j=1,nlobi)
359 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
360 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
367 c blower(k,l,j)=gaussc(ind,j,i)
372 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
373 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
380 C Read scrot parameters for potentials determined from all-atom AM1 calculations
381 C added by Urszula Kozlowska 07/11/2007
389 read(irotam,*) sc_parmin(j,i)
397 C Read torsional parameters in old format
399 read (itorp,*) ntortyp,nterm_old
400 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
401 read (itorp,*) (itortyp(i),i=1,ntyp)
406 read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
412 write (iout,'(/a/)') 'Torsional constants:'
415 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
416 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
422 C Read torsional parameters
424 read (itorp,*) ntortyp
425 read (itorp,*) (itortyp(i),i=1,ntyp)
426 write (iout,*) 'ntortyp',ntortyp
429 itortyp(i)=-itortyp(-i)
431 c write (iout,*) 'ntortyp',ntortyp
433 do j=-ntortyp+1,ntortyp-1
434 read (itorp,*) nterm(i,j,iblock),
436 nterm(-i,-j,iblock)=nterm(i,j,iblock)
437 nlor(-i,-j,iblock)=nlor(i,j,iblock)
440 do k=1,nterm(i,j,iblock)
441 read (itorp,*) kk,v1(k,i,j,iblock),
443 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
444 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
445 v0ij=v0ij+si*v1(k,i,j,iblock)
448 do k=1,nlor(i,j,iblock)
449 read (itorp,*) kk,vlor1(k,i,j),
450 & vlor2(k,i,j),vlor3(k,i,j)
451 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
454 v0(-i,-j,iblock)=v0ij
460 write (iout,'(/a/)') 'Torsional constants:'
463 write (iout,*) 'ityp',i,' jtyp',j
464 write (iout,*) 'Fourier constants'
465 do k=1,nterm(i,j,iblock)
466 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
469 write (iout,*) 'Lorenz constants'
470 do k=1,nlor(i,j,iblock)
471 write (iout,'(3(1pe15.5))')
472 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
478 C 6/23/01 Read parameters for double torsionals
482 do j=-ntortyp+1,ntortyp-1
483 do k=-ntortyp+1,ntortyp-1
484 read (itordp,'(3a1)') t1,t2,t3
485 if (t1.ne.onelett(i) .or. t2.ne.onelett(j)
486 & .or. t3.ne.onelett(k)) then
487 write (iout,*) "Error in double torsional parameter file",
489 stop "Error in double torsional parameter file"
491 read (itordp,*) ntermd_1(i,j,k,iblock),
492 & ntermd_2(i,j,k,iblock)
493 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
494 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
495 read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,
496 & ntermd_1(i,j,k,iblock))
497 read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,
498 & ntermd_1(i,j,k,iblock))
499 read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,
500 & ntermd_1(i,j,k,iblock))
501 read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,
502 & ntermd_1(i,j,k,iblock))
503 C Martix of D parameters for one dimesional foureir series
504 do l=1,ntermd_1(i,j,k,iblock)
505 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
506 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
507 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
508 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
509 c write(iout,*) "whcodze" ,
510 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
512 read (itordp,*) ((v2c(l,m,i,j,k,iblock),
513 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
514 & v2s(m,l,i,j,k,iblock),
515 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
516 C Martix of D parameters for two dimesional fourier series
517 do l=1,ntermd_2(i,j,k,iblock)
519 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
520 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
521 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
522 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
531 write (iout,*) 'Constants for double torsionals'
534 do j=-ntortyp+1,ntortyp-1
535 do k=-ntortyp+1,ntortyp-1
536 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
537 & ' nsingle',ntermd_1(i,j,k,iblock),
538 & ' ndouble',ntermd_2(i,j,k,iblock)
540 write (iout,*) 'Single angles:'
541 do l=1,ntermd_1(i,j,k,iblock)
542 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
543 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
544 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
545 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
548 write (iout,*) 'Pairs of angles:'
549 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
550 do l=1,ntermd_2(i,j,k,iblock)
551 write (iout,'(i5,20f10.5)')
552 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
555 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
556 do l=1,ntermd_2(i,j,k,iblock)
557 write (iout,'(i5,20f10.5)')
558 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
559 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
569 C Read of Side-chain backbone correlation parameters
570 C Modified 11 May 2012 by Adasko
572 read (isccor,*) nsccortyp
573 read (isccor,*) (isccortyp(i),i=1,ntyp)
575 isccortyp(i)=-isccortyp(-i)
577 iscprol=isccortyp(20)
578 c write (iout,*) 'ntortyp',ntortyp
580 cc maxinter is maximum interaction sites
585 &nterm_sccor(i,j),nlor_sccor(i,j)
591 nterm_sccor(-i,j)=nterm_sccor(i,j)
592 nterm_sccor(-i,-j)=nterm_sccor(i,j)
593 nterm_sccor(i,-j)=nterm_sccor(i,j)
594 do k=1,nterm_sccor(i,j)
595 read (isccor,*) kk,v1sccor(k,l,i,j)
597 if (j.eq.iscprol) then
598 if (i.eq.isccortyp(10)) then
599 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
600 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
602 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
603 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
604 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
605 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
606 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
607 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
608 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
609 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
612 if (i.eq.isccortyp(10)) then
613 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
614 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
616 if (j.eq.isccortyp(10)) then
617 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
618 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
620 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
621 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
622 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
623 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
624 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
625 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
629 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
630 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
631 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
632 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
635 do k=1,nlor_sccor(i,j)
636 read (isccor,*) kk,vlor1sccor(k,i,j),
637 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
638 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
639 &(1+vlor3sccor(k,i,j)**2)
641 v0sccor(l,i,j)=v0ijsccor
642 v0sccor(l,-i,j)=v0ijsccor1
643 v0sccor(l,i,-j)=v0ijsccor2
644 v0sccor(l,-i,-j)=v0ijsccor3
650 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
653 write (iout,*) 'ityp',i,' jtyp',j
654 write (iout,*) 'Fourier constants'
655 do k=1,nterm_sccor(i,j)
656 write (iout,'(2(1pe15.5))')
657 & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
659 write (iout,*) 'Lorenz constants'
660 do k=1,nlor_sccor(i,j)
661 write (iout,'(3(1pe15.5))')
662 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
668 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
669 C interaction energy of the Gly, Ala, and Pro prototypes.
671 read (ifourier,*) nloctyp
674 read (ifourier,*) (b(ii,i),ii=1,13)
676 write (iout,*) 'Type',i
677 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
683 B1tilde(1,i) = b(3,i)
684 B1tilde(2,i) =-b(5,i)
685 B1tilde(1,-i) =-b(3,i)
686 B1tilde(2,-i) =b(5,i)
701 Ctilde(2,1,i)=-b(9,i)
703 Ctilde(1,1,-i)=b(7,i)
704 Ctilde(1,2,-i)=-b(9,i)
705 Ctilde(2,1,-i)=b(9,i)
706 Ctilde(2,2,-i)=b(7,i)
717 Dtilde(2,1,i)=-b(8,i)
719 Dtilde(1,1,-i)=b(6,i)
720 Dtilde(1,2,-i)=-b(8,i)
721 Dtilde(2,1,-i)=b(8,i)
722 Dtilde(2,2,-i)=b(6,i)
723 EE(1,1,i)= b(10,i)+b(11,i)
724 EE(2,2,i)=-b(10,i)+b(11,i)
725 EE(2,1,i)= b(12,i)-b(13,i)
726 EE(1,2,i)= b(12,i)+b(13,i)
727 EE(1,1,-i)= b(10,i)+b(11,i)
728 EE(2,2,-i)=-b(10,i)+b(11,i)
729 EE(2,1,-i)=-b(12,i)+b(13,i)
730 EE(1,2,-i)=-b(12,i)-b(13,i)
734 write (iout,*) 'Type',i
736 c write (iout,'(f10.5)') B1(:,i)
737 write(iout,*) B1(1,i),B1(2,i)
739 c write (iout,'(f10.5)') B2(:,i)
740 write(iout,*) B2(1,i),B2(2,i)
743 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
747 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
751 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
756 C Read electrostatic-interaction parameters
759 write (iout,'(/a)') 'Electrostatic interaction constants:'
760 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
761 & 'IT','JT','APP','BPP','AEL6','AEL3'
763 read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
764 read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
765 read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
766 read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
771 app (i,j)=epp(i,j)*rri*rri
772 bpp (i,j)=-2.0D0*epp(i,j)*rri
773 ael6(i,j)=elpp6(i,j)*4.2D0**6
774 ael3(i,j)=elpp3(i,j)*4.2D0**3
775 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
776 & ael6(i,j),ael3(i,j)
780 C Read side-chain interaction parameters.
782 read (isidep,*) ipot,expon
783 if (ipot.lt.1 .or. ipot.gt.5) then
784 write (iout,'(2a)') 'Error while reading SC interaction',
785 & 'potential file - unknown potential type.'
789 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
790 & ', exponents are ',expon,2*expon
791 goto (10,20,30,30,40) ipot
792 C----------------------- LJ potential ---------------------------------
793 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
795 write (iout,'(/a/)') 'Parameters of the LJ potential:'
796 write (iout,'(a/)') 'The epsilon array:'
797 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
798 write (iout,'(/a)') 'One-body parameters:'
799 write (iout,'(a,4x,a)') 'residue','sigma'
800 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
803 C----------------------- LJK potential --------------------------------
804 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
805 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
807 write (iout,'(/a/)') 'Parameters of the LJK potential:'
808 write (iout,'(a/)') 'The epsilon array:'
809 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
810 write (iout,'(/a)') 'One-body parameters:'
811 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
812 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
816 C---------------------- GB or BP potential -----------------------------
817 30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
818 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
820 C For the GB potential convert sigma'**2 into chi'
823 chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
827 write (iout,'(/a/)') 'Parameters of the BP potential:'
828 write (iout,'(a/)') 'The epsilon array:'
829 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
830 write (iout,'(/a)') 'One-body parameters:'
831 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
833 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
834 & chip(i),alp(i),i=1,ntyp)
837 C--------------------- GBV potential -----------------------------------
838 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
839 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
840 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
842 write (iout,'(/a/)') 'Parameters of the GBV potential:'
843 write (iout,'(a/)') 'The epsilon array:'
844 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
845 write (iout,'(/a)') 'One-body parameters:'
846 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
847 & 's||/s_|_^2',' chip ',' alph '
848 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
849 & sigii(i),chip(i),alp(i),i=1,ntyp)
853 C-----------------------------------------------------------------------
854 C Calculate the "working" parameters of SC interactions.
862 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
863 sigma(j,i)=sigma(i,j)
864 rs0(i,j)=dwa16*sigma(i,j)
868 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
869 & 'Working parameters of the SC interactions:',
870 & ' a ',' b ',' augm ',' sigma ',' r0 ',
875 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
884 sigeps=dsign(1.0D0,epsij)
886 aa(i,j)=epsij*rrij*rrij
887 bb(i,j)=-sigeps*epsij*rrij
895 ratsig1=sigt2sq/sigt1sq
896 ratsig2=1.0D0/ratsig1
897 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
898 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
899 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
903 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
904 sigmaii(i,j)=rsum_max
905 sigmaii(j,i)=rsum_max
907 c sigmaii(i,j)=r0(i,j)
908 c sigmaii(j,i)=r0(i,j)
910 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
911 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
912 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
913 augm(i,j)=epsij*r_augm**(2*expon)
914 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
921 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
922 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
923 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
928 C Define the SC-p interaction constants
938 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
940 c aad(i,1)=0.3D0*4.0D0**12
941 C Following line for constants currently implemented
942 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
943 aad(i,1)=1.5D0*4.0D0**12
944 c aad(i,1)=0.17D0*5.6D0**12
946 C "Soft" SC-p repulsion
948 C Following line for constants currently implemented
949 c aad(i,1)=0.3D0*4.0D0**6
950 C "Hard" SC-p repulsion
951 bad(i,1)=3.0D0*4.0D0**6
952 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
961 C 8/9/01 Read the SC-p interaction constants from file
964 read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
967 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
968 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
969 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
970 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
974 write (iout,*) "Parameters of SC-p interactions:"
976 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
977 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
982 C Define the constants of the disulfide bridge
986 c Old arbitrary potential - commented out.
991 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
992 c energy surface of diethyl disulfide.
993 c A. Liwo and U. Kozlowska, 11/24/03
1003 write (iout,'(/a)') "Disulfide bridge parameters:"
1004 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1005 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1006 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1007 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,