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))
468 C Read of Side-chain backbone correlation parameters
469 C Modified 11 May 2012 by Adasko
472 read (isccor,*) nsccortyp
473 read (isccor,*) (isccortyp(i),i=1,ntyp)
474 c write (iout,*) 'ntortyp',ntortyp
476 cc maxinter is maximum interaction sites
481 &nterm_sccor(i,j),nlor_sccor(i,j)
485 do k=1,nterm_sccor(i,j)
486 read (isccor,*) kk,v1sccor(k,l,i,j)
488 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
491 do k=1,nlor_sccor(i,j)
492 read (isccor,*) kk,vlor1sccor(k,i,j),
493 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
494 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
495 &(1+vlor3sccor(k,i,j)**2)
497 v0sccor(i,j)=v0ijsccor
504 write (iout,'(/a/)') 'Torsional constants:'
508 write (iout,*) 'ityp',i,' jtyp',j
509 write (iout,*) 'Fourier constants'
510 do k=1,nterm_sccor(i,j)
511 write (iout,'(2(1pe15.5))')
512 & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
514 write (iout,*) 'Lorenz constants'
515 do k=1,nlor_sccor(i,j)
516 write (iout,'(3(1pe15.5))')
517 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
525 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
526 C interaction energy of the Gly, Ala, and Pro prototypes.
528 read (ifourier,*) nloctyp
531 read (ifourier,*) (b(ii,i),ii=1,13)
533 write (iout,*) 'Type',i
534 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
538 B1tilde(1,i) = b(3,i)
539 B1tilde(2,i) =-b(5,i)
548 Ctilde(2,1,i)=-b(9,i)
556 Dtilde(2,1,i)=-b(8,i)
558 EE(1,1,i)= b(10,i)+b(11,i)
559 EE(2,2,i)=-b(10,i)+b(11,i)
560 EE(2,1,i)= b(12,i)-b(13,i)
561 EE(1,2,i)= b(12,i)+b(13,i)
565 write (iout,*) 'Type',i
567 c write (iout,'(f10.5)') B1(:,i)
568 write(iout,*) B1(1,i),B1(2,i)
570 c write (iout,'(f10.5)') B2(:,i)
571 write(iout,*) B2(1,i),B2(2,i)
574 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
578 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
582 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
587 C Read electrostatic-interaction parameters
590 write (iout,'(/a)') 'Electrostatic interaction constants:'
591 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
592 & 'IT','JT','APP','BPP','AEL6','AEL3'
594 read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
595 read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
596 read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
597 read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
602 app (i,j)=epp(i,j)*rri*rri
603 bpp (i,j)=-2.0D0*epp(i,j)*rri
604 ael6(i,j)=elpp6(i,j)*4.2D0**6
605 ael3(i,j)=elpp3(i,j)*4.2D0**3
606 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
607 & ael6(i,j),ael3(i,j)
611 C Read side-chain interaction parameters.
613 read (isidep,*) ipot,expon
614 if (ipot.lt.1 .or. ipot.gt.5) then
615 write (iout,'(2a)') 'Error while reading SC interaction',
616 & 'potential file - unknown potential type.'
620 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
621 & ', exponents are ',expon,2*expon
622 goto (10,20,30,30,40) ipot
623 C----------------------- LJ potential ---------------------------------
624 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
626 write (iout,'(/a/)') 'Parameters of the LJ potential:'
627 write (iout,'(a/)') 'The epsilon array:'
628 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
629 write (iout,'(/a)') 'One-body parameters:'
630 write (iout,'(a,4x,a)') 'residue','sigma'
631 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
634 C----------------------- LJK potential --------------------------------
635 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
636 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
638 write (iout,'(/a/)') 'Parameters of the LJK potential:'
639 write (iout,'(a/)') 'The epsilon array:'
640 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
641 write (iout,'(/a)') 'One-body parameters:'
642 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
643 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
647 C---------------------- GB or BP potential -----------------------------
648 30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
649 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
651 C For the GB potential convert sigma'**2 into chi'
654 chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
658 write (iout,'(/a/)') 'Parameters of the BP potential:'
659 write (iout,'(a/)') 'The epsilon array:'
660 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
661 write (iout,'(/a)') 'One-body parameters:'
662 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
664 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
665 & chip(i),alp(i),i=1,ntyp)
668 C--------------------- GBV potential -----------------------------------
669 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
670 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
671 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
673 write (iout,'(/a/)') 'Parameters of the GBV potential:'
674 write (iout,'(a/)') 'The epsilon array:'
675 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
676 write (iout,'(/a)') 'One-body parameters:'
677 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
678 & 's||/s_|_^2',' chip ',' alph '
679 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
680 & sigii(i),chip(i),alp(i),i=1,ntyp)
684 C-----------------------------------------------------------------------
685 C Calculate the "working" parameters of SC interactions.
693 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
694 sigma(j,i)=sigma(i,j)
695 rs0(i,j)=dwa16*sigma(i,j)
699 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
700 & 'Working parameters of the SC interactions:',
701 & ' a ',' b ',' augm ',' sigma ',' r0 ',
706 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
715 sigeps=dsign(1.0D0,epsij)
717 aa(i,j)=epsij*rrij*rrij
718 bb(i,j)=-sigeps*epsij*rrij
726 ratsig1=sigt2sq/sigt1sq
727 ratsig2=1.0D0/ratsig1
728 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
729 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
730 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
734 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
735 sigmaii(i,j)=rsum_max
736 sigmaii(j,i)=rsum_max
738 c sigmaii(i,j)=r0(i,j)
739 c sigmaii(j,i)=r0(i,j)
741 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
742 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
743 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
744 augm(i,j)=epsij*r_augm**(2*expon)
745 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
752 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
753 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
754 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
759 C Define the SC-p interaction constants
769 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
771 c aad(i,1)=0.3D0*4.0D0**12
772 C Following line for constants currently implemented
773 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
774 aad(i,1)=1.5D0*4.0D0**12
775 c aad(i,1)=0.17D0*5.6D0**12
777 C "Soft" SC-p repulsion
779 C Following line for constants currently implemented
780 c aad(i,1)=0.3D0*4.0D0**6
781 C "Hard" SC-p repulsion
782 bad(i,1)=3.0D0*4.0D0**6
783 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
792 C 8/9/01 Read the SC-p interaction constants from file
795 read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
798 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
799 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
800 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
801 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
805 write (iout,*) "Parameters of SC-p interactions:"
807 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
808 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
813 C Define the constants of the disulfide bridge
817 c Old arbitrary potential - commented out.
822 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
823 c energy surface of diethyl disulfide.
824 c A. Liwo and U. Kozlowska, 11/24/03