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(4) /"G","A","P","D"/
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),j=1,2),(bthet(j,i),j=1,2)
82 read (ithep,*) (polthet(j,i),j=0,3)
83 read (ithep,*) (gthet(j,i),j=1,3)
84 read (ithep,*) theta0(i),sig0(i),sigc0(i)
90 c & 'Parameters of the virtual-bond valence angles:'
91 c write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
92 c & ' ATHETA0 ',' A1 ',' A2 ',
95 c write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
96 c & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
98 c write (iout,'(/a/9x,5a/79(1h-))')
99 c & 'Parameters of the expression for sigma(theta_c):',
100 c & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
101 c & ' ALPH3 ',' SIGMA0C '
103 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
104 c & (polthet(j,i),j=0,3),sigc0(i)
106 c write (iout,'(/a/9x,5a/79(1h-))')
107 c & 'Parameters of the second gaussian:',
108 c & ' THETA0 ',' SIGMA0 ',' G1 ',
111 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
112 c & sig0(i),(gthet(j,i),j=1,3)
115 & 'Parameters of the virtual-bond valence angles:'
116 write (iout,'(/a/9x,5a/79(1h-))')
117 & 'Coefficients of expansion',
118 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
119 & ' b1*10^1 ',' b2*10^1 '
121 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
122 & a0thet(i),(100*athet(j,i),j=1,2),(10*bthet(j,i),j=1,2)
124 write (iout,'(/a/9x,5a/79(1h-))')
125 & 'Parameters of the expression for sigma(theta_c):',
126 & ' alpha0 ',' alph1 ',' alph2 ',
127 & ' alhp3 ',' sigma0c '
129 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(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*10^2 ',' G1*10^-1',
137 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
138 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
143 C Read the parameters of Utheta determined from ab initio surfaces
144 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
146 read (ithep,*) nthetyp,ntheterm,ntheterm2,
147 & ntheterm3,nsingle,ndouble
148 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
149 read (ithep,*) (ithetyp(i),i=1,ntyp1)
155 aathet(l,i,j,k)=0.0d0
159 bbthet(m,l,i,j,k)=0.0d0
160 ccthet(m,l,i,j,k)=0.0d0
161 ddthet(m,l,i,j,k)=0.0d0
162 eethet(m,l,i,j,k)=0.0d0
168 ffthet(mm,m,l,i,j,k)=0.0d0
169 ggthet(mm,m,l,i,j,k)=0.0d0
179 read (ithep,'(3a)') res1,res2,res3
180 read (ithep,*) aa0thet(i,j,k)
181 read (ithep,*)(aathet(l,i,j,k),l=1,ntheterm)
183 & ((bbthet(lll,ll,i,j,k),lll=1,nsingle),
184 & (ccthet(lll,ll,i,j,k),lll=1,nsingle),
185 & (ddthet(lll,ll,i,j,k),lll=1,nsingle),
186 & (eethet(lll,ll,i,j,k),lll=1,nsingle),ll=1,ntheterm2)
188 & (((ffthet(llll,lll,ll,i,j,k),ffthet(lll,llll,ll,i,j,k),
189 & ggthet(llll,lll,ll,i,j,k),ggthet(lll,llll,ll,i,j,k),
190 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
195 C For dummy ends assign glycine-type coefficients of theta-only terms; the
196 C coefficients of theta-and-gamma-dependent terms are zero.
201 aathet(l,i,j,nthetyp+1)=aathet(l,i,j,1)
202 aathet(l,nthetyp+1,i,j)=aathet(l,1,i,j)
204 aa0thet(i,j,nthetyp+1)=aa0thet(i,j,1)
205 aa0thet(nthetyp+1,i,j)=aa0thet(1,i,j)
208 aathet(l,nthetyp+1,i,nthetyp+1)=aathet(l,1,i,1)
210 aa0thet(nthetyp+1,i,nthetyp+1)=aa0thet(1,i,1)
213 C Control printout of the coefficients of virtual-bond-angle potentials
216 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
220 write (iout,'(//4a)')
221 & 'Type ',onelett(i),onelett(j),onelett(k)
222 write (iout,'(//a,10x,a)') " l","a[l]"
223 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k)
224 write (iout,'(i2,1pe15.5)')
225 & (l,aathet(l,i,j,k),l=1,ntheterm)
227 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
228 & "b",l,"c",l,"d",l,"e",l
230 write (iout,'(i2,4(1pe15.5))') m,
231 & bbthet(m,l,i,j,k),ccthet(m,l,i,j,k),
232 & ddthet(m,l,i,j,k),eethet(m,l,i,j,k)
236 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
237 & "f+",l,"f-",l,"g+",l,"g-",l
240 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
241 & ffthet(n,m,l,i,j,k),ffthet(m,n,l,i,j,k),
242 & ggthet(n,m,l,i,j,k),ggthet(m,n,l,i,j,k)
255 C Read the parameters of the probability distribution/energy expression
256 C of the side chains.
259 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
263 dsc_inv(i)=1.0D0/dsc(i)
274 read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
276 read (irotam,*) bsc(j,i)
277 read (irotam,*) (censc(k,j,i),k=1,3),
278 & ((blower(k,l,j),l=1,k),k=1,3)
285 akl=akl+blower(k,m,j)*blower(l,m,j)
296 write (iout,'(/a)') 'Parameters of side-chain local geometry'
300 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
301 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
302 c write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
303 c write (iout,'(a,f10.4,4(16x,f10.4))')
304 c & 'Center ',(bsc(j,i),j=1,nlobi)
305 c write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
306 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
307 & 'log h',(bsc(j,i),j=1,nlobi)
308 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
309 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
316 c blower(k,l,j)=gaussc(ind,j,i)
321 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
322 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
329 C Read scrot parameters for potentials determined from all-atom AM1 calculations
330 C added by Urszula Kozlowska 07/11/2007
338 read(irotam,*) sc_parmin(j,i)
346 C Read torsional parameters in old format
348 read (itorp,*) ntortyp,nterm_old
349 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
350 read (itorp,*) (itortyp(i),i=1,ntyp)
355 read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
361 write (iout,'(/a/)') 'Torsional constants:'
364 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
365 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
371 C Read torsional parameters
373 read (itorp,*) ntortyp
374 read (itorp,*) (itortyp(i),i=1,ntyp)
375 write (iout,*) 'ntortyp',ntortyp
378 read (itorp,*) nterm(i,j),nlor(i,j)
382 read (itorp,*) kk,v1(k,i,j),v2(k,i,j)
383 v0ij=v0ij+si*v1(k,i,j)
387 read (itorp,*) kk,vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
388 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
395 write (iout,'(/a/)') 'Torsional constants:'
398 write (iout,*) 'ityp',i,' jtyp',j
399 write (iout,*) 'Fourier constants'
401 write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
403 write (iout,*) 'Lorenz constants'
405 write (iout,'(3(1pe15.5))')
406 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
412 C 6/23/01 Read parameters for double torsionals
417 read (itordp,'(3a1)') t1,t2,t3
418 if (t1.ne.onelett(i) .or. t2.ne.onelett(j)
419 & .or. t3.ne.onelett(k)) then
420 write (iout,*) "Error in double torsional parameter file",
422 stop "Error in double torsional parameter file"
424 read (itordp,*) ntermd_1(i,j,k),ntermd_2(i,j,k)
425 read (itordp,*) (v1c(1,l,i,j,k),l=1,ntermd_1(i,j,k))
426 read (itordp,*) (v1s(1,l,i,j,k),l=1,ntermd_1(i,j,k))
427 read (itordp,*) (v1c(2,l,i,j,k),l=1,ntermd_1(i,j,k))
428 read (itordp,*) (v1s(2,l,i,j,k),l=1,ntermd_1(i,j,k))
429 read (itordp,*) ((v2c(l,m,i,j,k),v2c(m,l,i,j,k),
430 & v2s(l,m,i,j,k),v2s(m,l,i,j,k),m=1,l-1),l=1,ntermd_2(i,j,k))
436 write (iout,*) 'Constants for double torsionals'
440 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
441 & ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k)
443 write (iout,*) 'Single angles:'
444 do l=1,ntermd_1(i,j,k)
445 write (iout,'(i5,2f10.5,5x,2f10.5)') l,
446 & v1c(1,l,i,j,k),v1s(1,l,i,j,k),
447 & v1c(2,l,i,j,k),v1s(2,l,i,j,k)
450 write (iout,*) 'Pairs of angles:'
451 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
452 do l=1,ntermd_2(i,j,k)
453 write (iout,'(i5,20f10.5)')
454 & l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k))
457 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
458 do l=1,ntermd_2(i,j,k)
459 write (iout,'(i5,20f10.5)')
460 & l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k))
469 C Read of Side-chain backbone correlation parameters
470 C Modified 11 May 2012 by Adasko
472 read (isccor,*) nsccortyp
473 read (isccor,*) (isccortyp(i),i=1,ntyp)
475 isccortyp(i)=-isccortyp(-i)
477 iscprol=isccortyp(20)
478 c write (iout,*) 'ntortyp',ntortyp
480 cc maxinter is maximum interaction sites
485 &nterm_sccor(i,j),nlor_sccor(i,j)
491 nterm_sccor(-i,j)=nterm_sccor(i,j)
492 nterm_sccor(-i,-j)=nterm_sccor(i,j)
493 nterm_sccor(i,-j)=nterm_sccor(i,j)
494 do k=1,nterm_sccor(i,j)
495 read (isccor,*) kk,v1sccor(k,l,i,j)
497 if (j.eq.iscprol) then
498 if (i.eq.isccortyp(10)) then
499 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
500 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
502 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
503 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
504 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
505 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
506 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
507 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
508 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
509 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
512 if (i.eq.isccortyp(10)) then
513 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
514 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
516 if (j.eq.isccortyp(10)) then
517 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
518 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
520 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
521 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
522 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
523 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
524 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
525 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
529 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
530 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
531 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
532 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
535 do k=1,nlor_sccor(i,j)
536 read (isccor,*) kk,vlor1sccor(k,i,j),
537 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
538 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
539 &(1+vlor3sccor(k,i,j)**2)
541 v0sccor(l,i,j)=v0ijsccor
542 v0sccor(l,-i,j)=v0ijsccor1
543 v0sccor(l,i,-j)=v0ijsccor2
544 v0sccor(l,-i,-j)=v0ijsccor3
550 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
554 write (iout,*) 'ityp',i,' jtyp',j
555 write (iout,*) 'Fourier constants'
556 do k=1,nterm_sccor(i,j)
557 write (iout,'(2(1pe15.5))')
558 & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
560 write (iout,*) 'Lorenz constants'
561 do k=1,nlor_sccor(i,j)
562 write (iout,'(3(1pe15.5))')
563 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
570 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
571 C interaction energy of the Gly, Ala, and Pro prototypes.
573 read (ifourier,*) nloctyp
576 read (ifourier,*) (b(ii,i),ii=1,13)
578 write (iout,*) 'Type',i
579 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
583 B1tilde(1,i) = b(3,i)
584 B1tilde(2,i) =-b(5,i)
593 Ctilde(2,1,i)=-b(9,i)
601 Dtilde(2,1,i)=-b(8,i)
603 EE(1,1,i)= b(10,i)+b(11,i)
604 EE(2,2,i)=-b(10,i)+b(11,i)
605 EE(2,1,i)= b(12,i)-b(13,i)
606 EE(1,2,i)= b(12,i)+b(13,i)
610 write (iout,*) 'Type',i
612 c write (iout,'(f10.5)') B1(:,i)
613 write(iout,*) B1(1,i),B1(2,i)
615 c write (iout,'(f10.5)') B2(:,i)
616 write(iout,*) B2(1,i),B2(2,i)
619 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
623 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
627 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
632 C Read electrostatic-interaction parameters
635 write (iout,'(/a)') 'Electrostatic interaction constants:'
636 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
637 & 'IT','JT','APP','BPP','AEL6','AEL3'
639 read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
640 read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
641 read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
642 read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
647 app (i,j)=epp(i,j)*rri*rri
648 bpp (i,j)=-2.0D0*epp(i,j)*rri
649 ael6(i,j)=elpp6(i,j)*4.2D0**6
650 ael3(i,j)=elpp3(i,j)*4.2D0**3
651 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
652 & ael6(i,j),ael3(i,j)
656 C Read side-chain interaction parameters.
658 read (isidep,*) ipot,expon
659 if (ipot.lt.1 .or. ipot.gt.5) then
660 write (iout,'(2a)') 'Error while reading SC interaction',
661 & 'potential file - unknown potential type.'
665 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
666 & ', exponents are ',expon,2*expon
667 goto (10,20,30,30,40) ipot
668 C----------------------- LJ potential ---------------------------------
669 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
671 write (iout,'(/a/)') 'Parameters of the LJ potential:'
672 write (iout,'(a/)') 'The epsilon array:'
673 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
674 write (iout,'(/a)') 'One-body parameters:'
675 write (iout,'(a,4x,a)') 'residue','sigma'
676 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
679 C----------------------- LJK potential --------------------------------
680 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
681 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
683 write (iout,'(/a/)') 'Parameters of the LJK potential:'
684 write (iout,'(a/)') 'The epsilon array:'
685 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
686 write (iout,'(/a)') 'One-body parameters:'
687 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
688 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
692 C---------------------- GB or BP potential -----------------------------
693 30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
694 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
696 C For the GB potential convert sigma'**2 into chi'
699 chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
703 write (iout,'(/a/)') 'Parameters of the BP potential:'
704 write (iout,'(a/)') 'The epsilon array:'
705 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
706 write (iout,'(/a)') 'One-body parameters:'
707 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
709 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
710 & chip(i),alp(i),i=1,ntyp)
713 C--------------------- GBV potential -----------------------------------
714 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
715 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
716 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
718 write (iout,'(/a/)') 'Parameters of the GBV potential:'
719 write (iout,'(a/)') 'The epsilon array:'
720 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
721 write (iout,'(/a)') 'One-body parameters:'
722 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
723 & 's||/s_|_^2',' chip ',' alph '
724 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
725 & sigii(i),chip(i),alp(i),i=1,ntyp)
729 C-----------------------------------------------------------------------
730 C Calculate the "working" parameters of SC interactions.
738 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
739 sigma(j,i)=sigma(i,j)
740 rs0(i,j)=dwa16*sigma(i,j)
744 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
745 & 'Working parameters of the SC interactions:',
746 & ' a ',' b ',' augm ',' sigma ',' r0 ',
751 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
760 sigeps=dsign(1.0D0,epsij)
762 aa(i,j)=epsij*rrij*rrij
763 bb(i,j)=-sigeps*epsij*rrij
771 ratsig1=sigt2sq/sigt1sq
772 ratsig2=1.0D0/ratsig1
773 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
774 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
775 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
779 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
780 sigmaii(i,j)=rsum_max
781 sigmaii(j,i)=rsum_max
783 c sigmaii(i,j)=r0(i,j)
784 c sigmaii(j,i)=r0(i,j)
786 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
787 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
788 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
789 augm(i,j)=epsij*r_augm**(2*expon)
790 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
797 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
798 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
799 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
804 C Define the SC-p interaction constants
814 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
816 c aad(i,1)=0.3D0*4.0D0**12
817 C Following line for constants currently implemented
818 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
819 aad(i,1)=1.5D0*4.0D0**12
820 c aad(i,1)=0.17D0*5.6D0**12
822 C "Soft" SC-p repulsion
824 C Following line for constants currently implemented
825 c aad(i,1)=0.3D0*4.0D0**6
826 C "Hard" SC-p repulsion
827 bad(i,1)=3.0D0*4.0D0**6
828 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
837 C 8/9/01 Read the SC-p interaction constants from file
840 read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
843 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
844 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
845 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
846 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
850 write (iout,*) "Parameters of SC-p interactions:"
852 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
853 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
858 C Define the constants of the disulfide bridge
862 c Old arbitrary potential - commented out.
867 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
868 c energy surface of diethyl disulfide.
869 c A. Liwo and U. Kozlowska, 11/24/03
879 write (iout,'(/a)') "Disulfide bridge parameters:"
880 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
881 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
882 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
883 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,