1 subroutine parmread(iparm,*)
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)
8 include 'DIMENSIONS.ZSCOPT'
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'
20 include 'COMMON.WEIGHTS'
21 include 'COMMON.OPTIM'
22 include 'COMMON.ENERGIES'
24 character*1 onelett(4) /"G","A","P","D"/
26 dimension blower(3,3,maxlob)
27 double precision ip,mp
28 character*800 controlcard
29 character*256 bondname_t,thetname_t,rotname_t,torname_t,
30 & tordname_t,fouriername_t,elename_t,sidename_t,scpname_t,
38 C Set LPRINT=.TRUE. for debugging
39 dwa16=2.0d0**(1.0d0/6.0d0)
42 C Read parameter filenames and open the parameter files.
44 write (iout,*) "iparm",iparm," myparm",myparm
46 call card_concat(controlcard,.true.)
48 if (iparm.eq.myparm) then
52 key = wname(i)(:ilen(wname(i)))
53 call reada(controlcard,key(:ilen(key)),ww(i),ww0(i))
58 c If reading not own parameters, skip assignment
60 call card_concat(controlcard,.false.)
62 if (iparm.ne.myparm) return
64 call reads(controlcard,"BONDPAR",bondname_t,bondname)
65 open (ibond,file=bondname_t,status='old')
67 call reads(controlcard,"THETPAR",thetname_t,thetname)
68 open (ithep,file=thetname_t,status='old')
70 call reads(controlcard,"ROTPAR",rotname_t,rotname)
71 open (irotam,file=rotname_t,status='old')
73 call reads(controlcard,"TORPAR",torname_t,torname)
74 open (itorp,file=torname_t,status='old')
76 call reads(controlcard,"TORDPAR",tordname_t,tordname)
77 open (itordp,file=tordname_t,status='old')
79 call reads(controlcard,"SCCORAR",sccorname_t,sccorname)
80 open (isccor,file=sccorname_t,status='old')
82 call reads(controlcard,"FOURIER",fouriername_t,fouriername)
83 open (ifourier,file=fouriername_t,status='old')
85 call reads(controlcard,"ELEPAR",elename_t,elename)
86 open (ielep,file=elename_t,status='old')
88 call reads(controlcard,"SIDEPAR",sidename_t,sidename)
89 open (isidep,file=sidename_t,status='old')
91 call reads(controlcard,"SCPPAR",scpname_t,scpname)
92 open (iscpp,file=scpname_t,status='old')
94 write (iout,*) "Parameter set:",iparm
95 write (iout,*) "Energy-term weights:"
96 write (iout,'(/5a12)') "WEIGHT ","START VALUE ",
97 & " MASK ","LOWER LIMIT ","UPPER LIMIT "
99 write (iout,'(a12,f12.5,i7,5x,2(f11.5,1x))') wname(i),ww(i),
100 & imask(i),ww_low(i),ww_up(i)
102 write (iout,*) "Sidechain potential file : ",
103 & sidename_t(:ilen(sidename_t))
105 write (iout,*) "SCp potential file : ",
106 & scpname_t(:ilen(scpname_t))
108 write (iout,*) "Electrostatic potential file : ",
109 & elename_t(:ilen(elename_t))
110 write (iout,*) "Cumulant coefficient file : ",
111 & fouriername_t(:ilen(fouriername_t))
112 write (iout,*) "Torsional parameter file : ",
113 & torname_t(:ilen(torname_t))
114 write (iout,*) "Double torsional parameter file : ",
115 & tordname_t(:ilen(tordname_t))
116 write (iout,*) "Backbone-rotamer parameter file : ",
117 & sccorname(:ilen(sccorname))
118 write (iout,*) "Bond & inertia constant file : ",
119 & bondname_t(:ilen(bondname_t))
120 write (iout,*) "Bending parameter file : ",
121 & thetname_t(:ilen(thetname_t))
122 write (iout,*) "Rotamer parameter file : ",
123 & rotname_t(:ilen(rotname_t))
125 C Assign virtual-bond length
128 vblinv2=vblinv*vblinv
130 read (ibond,*) vbldp0,akp
133 read (ibond,*) vbldsc0(1,i),aksc(1,i)
134 dsc(i) = vbldsc0(1,i)
138 dsc_inv(i)=1.0D0/dsc(i)
142 read (ibond,*) ijunk,vbldp0,akp,rjunk
144 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
146 dsc(i) = vbldsc0(1,i)
150 dsc_inv(i)=1.0D0/dsc(i)
155 write(iout,'(/a/)')"Force constants virtual bonds:"
156 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
158 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
160 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
161 & vbldsc0(1,i),aksc(1,i),abond0(1,i)
163 write (iout,'(13x,3f10.5)')
164 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
170 C Read the parameters of the probability distribution/energy expression
171 C of the virtual-bond valence angles theta
174 read (ithep,*) a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
175 read (ithep,*) (polthet(j,i),j=0,3)
176 read (ithep,*) (gthet(j,i),j=1,3)
177 read (ithep,*) theta0(i),sig0(i),sigc0(i)
183 c & 'Parameters of the virtual-bond valence angles:'
184 c write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
185 c & ' ATHETA0 ',' A1 ',' A2 ',
188 c write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
189 c & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
191 c write (iout,'(/a/9x,5a/79(1h-))')
192 c & 'Parameters of the expression for sigma(theta_c):',
193 c & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
194 c & ' ALPH3 ',' SIGMA0C '
196 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
197 c & (polthet(j,i),j=0,3),sigc0(i)
199 c write (iout,'(/a/9x,5a/79(1h-))')
200 c & 'Parameters of the second gaussian:',
201 c & ' THETA0 ',' SIGMA0 ',' G1 ',
204 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
205 c & sig0(i),(gthet(j,i),j=1,3)
208 & 'Parameters of the virtual-bond valence angles:'
209 write (iout,'(/a/9x,5a/79(1h-))')
210 & 'Coefficients of expansion',
211 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
212 & ' b1*10^1 ',' b2*10^1 '
214 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
215 & a0thet(i),(100*athet(j,i),j=1,2),(10*bthet(j,i),j=1,2)
217 write (iout,'(/a/9x,5a/79(1h-))')
218 & 'Parameters of the expression for sigma(theta_c):',
219 & ' alpha0 ',' alph1 ',' alph2 ',
220 & ' alhp3 ',' sigma0c '
222 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
223 & (polthet(j,i),j=0,3),sigc0(i)
225 write (iout,'(/a/9x,5a/79(1h-))')
226 & 'Parameters of the second gaussian:',
227 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
230 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
231 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
236 C Read the parameters of Utheta determined from ab initio surfaces
237 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
239 read (ithep,*) nthetyp,ntheterm,ntheterm2,
240 & ntheterm3,nsingle,ndouble
241 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
242 read (ithep,*) (ithetyp(i),i=1,ntyp1)
248 aathet(l,i,j,k)=0.0d0
252 bbthet(m,l,i,j,k)=0.0d0
253 ccthet(m,l,i,j,k)=0.0d0
254 ddthet(m,l,i,j,k)=0.0d0
255 eethet(m,l,i,j,k)=0.0d0
261 ffthet(mm,m,l,i,j,k)=0.0d0
262 ggthet(mm,m,l,i,j,k)=0.0d0
272 read (ithep,'(3a)') res1,res2,res3
273 read (ithep,*) aa0thet(i,j,k)
274 read (ithep,*)(aathet(l,i,j,k),l=1,ntheterm)
276 & ((bbthet(lll,ll,i,j,k),lll=1,nsingle),
277 & (ccthet(lll,ll,i,j,k),lll=1,nsingle),
278 & (ddthet(lll,ll,i,j,k),lll=1,nsingle),
279 & (eethet(lll,ll,i,j,k),lll=1,nsingle),ll=1,ntheterm2)
281 & (((ffthet(llll,lll,ll,i,j,k),ffthet(lll,llll,ll,i,j,k),
282 & ggthet(llll,lll,ll,i,j,k),ggthet(lll,llll,ll,i,j,k),
283 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
288 C For dummy ends assign glycine-type coefficients of theta-only terms; the
289 C coefficients of theta-and-gamma-dependent terms are zero.
294 aathet(l,i,j,nthetyp+1)=aathet(l,i,j,1)
295 aathet(l,nthetyp+1,i,j)=aathet(l,1,i,j)
297 aa0thet(i,j,nthetyp+1)=aa0thet(i,j,1)
298 aa0thet(nthetyp+1,i,j)=aa0thet(1,i,j)
301 aathet(l,nthetyp+1,i,nthetyp+1)=aathet(l,1,i,1)
303 aa0thet(nthetyp+1,i,nthetyp+1)=aa0thet(1,i,1)
306 C Control printout of the coefficients of virtual-bond-angle potentials
309 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
313 write (iout,'(//4a)')
314 & 'Type ',onelett(i),onelett(j),onelett(k)
315 write (iout,'(//a,10x,a)') " l","a[l]"
316 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k)
317 write (iout,'(i2,1pe15.5)')
318 & (l,aathet(l,i,j,k),l=1,ntheterm)
320 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
321 & "b",l,"c",l,"d",l,"e",l
323 write (iout,'(i2,4(1pe15.5))') m,
324 & bbthet(m,l,i,j,k),ccthet(m,l,i,j,k),
325 & ddthet(m,l,i,j,k),eethet(m,l,i,j,k)
329 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
330 & "f+",l,"f-",l,"g+",l,"g-",l
333 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
334 & ffthet(n,m,l,i,j,k),ffthet(m,n,l,i,j,k),
335 & ggthet(n,m,l,i,j,k),ggthet(m,n,l,i,j,k)
348 C Read the parameters of the probability distribution/energy expression
349 C of the side chains.
352 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
356 dsc_inv(i)=1.0D0/dsc(i)
367 read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
369 read (irotam,*) bsc(j,i)
370 read (irotam,*) (censc(k,j,i),k=1,3),
371 & ((blower(k,l,j),l=1,k),k=1,3)
378 akl=akl+blower(k,m,j)*blower(l,m,j)
389 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 c write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
396 c write (iout,'(a,f10.4,4(16x,f10.4))')
397 c & 'Center ',(bsc(j,i),j=1,nlobi)
398 c write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
399 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
400 & 'log h',(bsc(j,i),j=1,nlobi)
401 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
402 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
409 c blower(k,l,j)=gaussc(ind,j,i)
414 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
415 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
422 C Read scrot parameters for potentials determined from all-atom AM1 calculations
423 C added by Urszula Kozlowska 07/11/2007
431 read(irotam,*) sc_parmin(j,i)
439 C Read torsional parameters in old format
441 read (itorp,*) ntortyp,nterm_old
442 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
443 read (itorp,*) (itortyp(i),i=1,ntyp)
448 read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
454 write (iout,'(/a/)') 'Torsional constants:'
457 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
458 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
464 C Read torsional parameters
466 read (itorp,*) ntortyp
467 read (itorp,*) (itortyp(i),i=1,ntyp)
468 write (iout,*) 'ntortyp',ntortyp
471 read (itorp,*) nterm(i,j),nlor(i,j)
475 read (itorp,*) kk,v1(k,i,j),v2(k,i,j)
476 v0ij=v0ij+si*v1(k,i,j)
480 read (itorp,*) kk,vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
481 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
488 write (iout,'(/a/)') 'Torsional constants:'
491 write (iout,*) 'ityp',i,' jtyp',j
492 write (iout,*) 'Fourier constants'
494 write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
496 write (iout,*) 'Lorenz constants'
498 write (iout,'(3(1pe15.5))')
499 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
505 C 6/23/01 Read parameters for double torsionals
510 read (itordp,'(3a1)') t1,t2,t3
511 if (t1.ne.onelett(i) .or. t2.ne.onelett(j)
512 & .or. t3.ne.onelett(k)) then
513 write (iout,*) "Error in double torsional parameter file",
515 stop "Error in double torsional parameter file"
517 read (itordp,*) ntermd_1(i,j,k),ntermd_2(i,j,k)
518 read (itordp,*) (v1c(1,l,i,j,k),l=1,ntermd_1(i,j,k))
519 read (itordp,*) (v1s(1,l,i,j,k),l=1,ntermd_1(i,j,k))
520 read (itordp,*) (v1c(2,l,i,j,k),l=1,ntermd_1(i,j,k))
521 read (itordp,*) (v1s(2,l,i,j,k),l=1,ntermd_1(i,j,k))
522 read (itordp,*) ((v2c(l,m,i,j,k),v2c(m,l,i,j,k),
523 & v2s(l,m,i,j,k),v2s(m,l,i,j,k),m=1,l-1),l=1,ntermd_2(i,j,k))
529 write (iout,*) 'Constants for double torsionals'
533 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
534 & ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k)
536 write (iout,*) 'Single angles:'
537 do l=1,ntermd_1(i,j,k)
538 write (iout,'(i5,2f10.5,5x,2f10.5)') l,
539 & v1c(1,l,i,j,k),v1s(1,l,i,j,k),
540 & v1c(2,l,i,j,k),v1s(2,l,i,j,k)
543 write (iout,*) 'Pairs of angles:'
544 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
545 do l=1,ntermd_2(i,j,k)
546 write (iout,'(i5,20f10.5)')
547 & l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k))
550 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
551 do l=1,ntermd_2(i,j,k)
552 write (iout,'(i5,20f10.5)')
553 & l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k))
561 C Read of Side-chain backbone correlation parameters
562 C Modified 11 May 2012 by Adasko
565 read (isccor,*,end=1113,err=1113) nsccortyp
566 read (isccor,*,end=1113,err=1113) (isccortyp(i),i=1,ntyp)
567 write (iout,*) 'nsccortyp',nsccortyp
569 cc maxinter is maximum interaction sites
573 read (isccor,*,end=1113,err=1113) nterm_sccor(i,j),
578 DO k = 1, nterm_sccor(i,j)
579 read (isccor,*,end=1113,err=1113) kk, v1sccor(k,l,i,j)
581 v0ijsccor = v0ijsccor + si * v1sccor(k,l,i,j)
584 DO k = 1, nlor_sccor(i,j)
585 read (isccor,*,end=1113,err=1113) kk, vlor1sccor(k,i,j),
586 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
587 v0ijsccor = v0ijsccor + vlor1sccor(k,i,j) /
588 & (1 + vlor3sccor(k,i,j)**2)
590 v0sccor(i,j) = v0ijsccor
597 write (iout,'(/a/)') 'Torsional constants:'
601 write (iout,*) 'ityp',i,' jtyp',j
602 write (iout,*) 'Fourier constants'
603 DO k=1,nterm_sccor(i,j)
604 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
606 write (iout,*) 'Lorenz constants'
607 DO k=1,nlor_sccor(i,j)
608 write (iout,'(3(1pe15.5))')
609 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
616 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
617 C interaction energy of the Gly, Ala, and Pro prototypes.
619 read (ifourier,*) nloctyp
622 read (ifourier,*) (b(ii,i),ii=1,13)
624 write (iout,*) 'Type',i
625 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
629 B1tilde(1,i) = b(3,i)
630 B1tilde(2,i) =-b(5,i)
639 Ctilde(2,1,i)=-b(9,i)
647 Dtilde(2,1,i)=-b(8,i)
649 EE(1,1,i)= b(10,i)+b(11,i)
650 EE(2,2,i)=-b(10,i)+b(11,i)
651 EE(2,1,i)= b(12,i)-b(13,i)
652 EE(1,2,i)= b(12,i)+b(13,i)
656 write (iout,*) 'Type',i
658 c write (iout,'(f10.5)') B1(:,i)
659 write(iout,*) B1(1,i),B1(2,i)
661 c write (iout,'(f10.5)') B2(:,i)
662 write(iout,*) B2(1,i),B2(2,i)
665 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
669 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
673 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
678 C Read electrostatic-interaction parameters
681 write (iout,'(/a)') 'Electrostatic interaction constants:'
682 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
683 & 'IT','JT','APP','BPP','AEL6','AEL3'
685 read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
686 read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
687 read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
688 read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
693 app (i,j)=epp(i,j)*rri*rri
694 bpp (i,j)=-2.0D0*epp(i,j)*rri
695 ael6(i,j)=elpp6(i,j)*4.2D0**6
696 ael3(i,j)=elpp3(i,j)*4.2D0**3
697 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
698 & ael6(i,j),ael3(i,j)
702 C Read side-chain interaction parameters.
704 read (isidep,*) ipot,expon
705 if (ipot.lt.1 .or. ipot.gt.6) then
706 write (iout,'(2a)') 'Error while reading SC interaction',
707 & 'potential file - unknown potential type.'
711 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
712 & ', exponents are ',expon,2*expon
713 goto (10,20,30,30,40,50,60) ipot
714 C----------------------- LJ potential ---------------------------------
715 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
717 write (iout,'(/a/)') 'Parameters of the LJ potential:'
718 write (iout,'(a/)') 'The epsilon array:'
719 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
720 write (iout,'(/a)') 'One-body parameters:'
721 write (iout,'(a,4x,a)') 'residue','sigma'
722 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
725 C----------------------- LJK potential --------------------------------
726 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
727 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
729 write (iout,'(/a/)') 'Parameters of the LJK potential:'
730 write (iout,'(a/)') 'The epsilon array:'
731 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
732 write (iout,'(/a)') 'One-body parameters:'
733 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
734 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
738 C---------------------- GB or BP potential -----------------------------
739 30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
740 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
742 C For the GB potential convert sigma'**2 into chi'
745 chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
749 write (iout,'(/a/)') 'Parameters of the BP potential:'
750 write (iout,'(a/)') 'The epsilon array:'
751 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
752 write (iout,'(/a)') 'One-body parameters:'
753 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
755 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
756 & chip(i),alp(i),i=1,ntyp)
759 C--------------------- GBV potential -----------------------------------
760 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
761 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
762 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
764 write (iout,'(/a/)') 'Parameters of the GBV potential:'
765 write (iout,'(a/)') 'The epsilon array:'
766 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
767 write (iout,'(/a)') 'One-body parameters:'
768 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
769 & 's||/s_|_^2',' chip ',' alph '
770 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
771 & sigii(i),chip(i),alp(i),i=1,ntyp)
774 C--------------------- Momo potential -----------------------------------
776 read (isidep,*) (icharge(i),i=1,ntyp)
777 c write (2,*) "icharge",(icharge(i),i=1,ntyp)
780 c! write (*,*) "Im in ", i, " ", j
782 & eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i),
783 & (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j),
784 & chis(i,j),chis(j,i),
785 & nstate(i,j),(wstate(k,i,j),k=1,4),
790 & dtail(1,i,j),dtail(2,i,j),
791 & epshead(i,j),sig0head(i,j),
792 & rborn(i,j),rborn(j,i),
793 & (wqdip(k,i,j),k=1,2),wquad(i,j),
794 & alphapol(i,j),alphapol(j,i),
795 & (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j)
798 c! write (*,*) "nstate gly-gly", nstate(10,10)
799 c! THIS LOOP FILLS PARAMETERS FOR PAIRS OF AA's NOT EXPLICITLY
800 c! DEFINED IN SCPARM.MOMO. IT DOES SO BY TAKING THEM FROM SYMMETRIC
801 c! PAIR, FOR EG. IF ARG-HIS IS BLANK, IT WILL TAKE PARAMETERS
804 c! DISABLE IT AT >>YOUR OWN PERIL<<
809 sigma(i,j) = sigma(j,i)
810 nstate(i,j) = nstate(j,i)
811 sigmap1(i,j) = sigmap1(j,i)
812 sigmap2(i,j) = sigmap2(j,i)
813 sigiso1(i,j) = sigiso1(j,i)
814 sigiso2(i,j) = sigiso2(j,i)
817 alphasur(k,i,j) = alphasur(k,j,i)
818 wstate(k,i,j) = wstate(k,j,i)
819 alphiso(k,i,j) = alphiso(k,j,i)
822 dhead(2,1,i,j) = dhead(1,1,j,i)
823 dhead(2,2,i,j) = dhead(1,2,j,i)
824 dhead(1,1,i,j) = dhead(2,1,j,i)
825 dhead(1,2,i,j) = dhead(2,2,j,i)
826 dtail(1,i,j) = dtail(1,j,i)
827 dtail(2,i,j) = dtail(2,j,i)
830 c! write (*,*) "dhead(k,l,j,i) = ", dhead(k,l,j,i)
831 c! write (*,*) "dhead(k,l,i,j) = ", dhead(k,l,i,j)
832 c! dhead(l,k,i,j) = dhead(k,l,j,i)
836 epshead(i,j) = epshead(j,i)
837 sig0head(i,j) = sig0head(j,i)
840 wqdip(k,i,j) = wqdip(k,j,i)
843 wquad(i,j) = wquad(j,i)
844 epsintab(i,j) = epsintab(j,i)
849 if (.not.lprint) goto 70
851 & "Parameters of the new physics-based SC-SC interaction potential"
852 write (iout,'(/7a)') 'Residues',' epsGB',' rGB',
853 & ' chi1GB',' chi2GB',' chip1GB',' chip2GB'
856 write (iout,'(2(a3,1x),1pe10.3,5(0pf10.3))')
857 & restyp(i),restyp(j),eps(i,j),sigma(i,j),chi(i,j),chi(j,i),
858 & chipp(i,j),chipp(j,i)
861 write (iout,'(/9a)') 'Residues',' alphasur1',' alphasur2',
862 & ' alphasur3',' alphasur4',' sigmap1',' sigmap2',
866 write (iout,'(2(a3,1x),8(0pf10.3))')
867 & restyp(i),restyp(j),(alphasur(k,i,j),k=1,4),
868 & sigmap1(i,j),sigmap2(j,i),chis(i,j),chis(j,i)
871 write (iout,'(/14a)') 'Residues',' nst ',' wst1',
872 & ' wst2',' wst3',' wst4',' dh11',' dh21',
873 & ' dh12',' dh22',' dt1',' dt2',' epsh1',
877 write (iout,'(2(a3,1x),i3,4f8.3,6f7.2,f9.5,f7.2)')
878 & restyp(i),restyp(j),nstate(i,j),(wstate(k,i,j),k=1,4),
879 & ((dhead(l,k,i,j),l=1,2),k=1,2),dtail(1,i,j),dtail(2,i,j),
880 & epshead(i,j),sig0head(i,j)
883 write (iout,'(/12a)') 'Residues',' ch1',' ch2',
884 & ' rborn1',' rborn2',' wqdip1',' wqdip2',
888 write (iout,'(2(a3,1x),2i4,5f10.3)')
889 & restyp(i),restyp(j),icharge(i),icharge(j),
890 & rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j)
893 write (iout,'(/12a)') 'Residues',
895 & ' alphpol2',' alphiso1',' alpiso2',
896 & ' alpiso3',' alpiso4',' sigiso1',' sigiso2',
900 write (iout,'(2(a3,1x),11f10.3)')
901 & restyp(i),restyp(j),alphapol(i,j),alphapol(j,i),
902 & (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(j,i),
907 C For the GB potential convert sigma'**2 into chi'
909 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
912 write (iout,'(/a/)') 'Parameters of the BP potential:'
913 write (iout,'(a/)') 'The epsilon array:'
914 CALL printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
915 write (iout,'(/a)') 'One-body parameters:'
916 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
918 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
919 & chip(i),alp(i),i=1,ntyp)
923 C-----------------------------------------------------------------------
924 C Calculate the "working" parameters of SC interactions.
929 alphasur(k,j,i)=alphasur(k,i,j)
930 alphiso(k,j,i)=alphiso(k,i,j)
931 wstate(k,j,i)=wstate(k,i,j)
934 wqdip(k,j,i)=wqdip(k,i,j)
938 dhead(l,k,j,i)=dhead(l,k,i,j)
946 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
947 sigma(j,i)=sigma(i,j)
948 rs0(i,j)=dwa16*sigma(i,j)
955 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
956 & 'Working parameters of the SC interactions:',
957 & ' a ',' b ',' augm ',' sigma ',' r0 ',
962 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4 .or. ipot.eq.6 ) THEN
971 sigeps=dsign(1.0D0,epsij)
973 aa(i,j)=epsij*rrij*rrij
974 bb(i,j)=-sigeps*epsij*rrij
977 IF ((ipot.gt.2).AND.(ipot.LT.6)) THEN
982 ratsig1=sigt2sq/sigt1sq
983 ratsig2=1.0D0/ratsig1
984 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
985 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
986 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
990 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
991 sigmaii(i,j)=rsum_max
992 sigmaii(j,i)=rsum_max
994 c sigmaii(i,j)=r0(i,j)
995 c sigmaii(j,i)=r0(i,j)
997 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
998 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
999 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1000 augm(i,j)=epsij*r_augm**(2*expon)
1001 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1009 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1010 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1011 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1013 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3),2i3,10f8.4,
1015 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1016 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i),
1017 & icharge(i),icharge(j),chipp(i,j),chipp(j,i),
1018 & (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(j,i),
1019 & chis(i,j),chis(j,i),
1020 & nstate(i,j),(wstate(k,i,j),k=1,4),
1021 & ((dhead(l,k,i,j),l=1,2),k=1,2),dtail(1,i,j),dtail(2,i,j),
1022 & epshead(i,j),sig0head(i,j),
1023 & rborn(i,j),(wqdip(k,i,j),k=1,2),wquad(i,j),
1024 & alphapol(i,j),alphapol(j,i),
1025 & (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j)
1033 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1035 c aad(i,1)=0.3D0*4.0D0**12
1036 C Following line for constants currently implemented
1037 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1038 aad(i,1)=1.5D0*4.0D0**12
1039 c aad(i,1)=0.17D0*5.6D0**12
1041 C "Soft" SC-p repulsion
1043 C Following line for constants currently implemented
1044 c aad(i,1)=0.3D0*4.0D0**6
1045 C "Hard" SC-p repulsion
1046 bad(i,1)=3.0D0*4.0D0**6
1047 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1056 C 8/9/01 Read the SC-p interaction constants from file
1059 read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1062 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1063 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1064 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1065 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1069 write (iout,*) "Parameters of SC-p interactions:"
1071 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1072 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1077 C Define the constants of the disulfide bridge
1081 c Old arbitrary potential - commented out.
1086 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1087 c energy surface of diethyl disulfide.
1088 c A. Liwo and U. Kozlowska, 11/24/03
1098 write (iout,'(/a)') "Disulfide bridge parameters:"
1099 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1100 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1101 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1102 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1106 & "Error reading side-chain torsional energy parameters."