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 'DIMENSIONS.FREE'
10 include 'COMMON.IOUNITS'
11 include 'COMMON.CHAIN'
12 include 'COMMON.INTERACT'
14 include 'COMMON.LOCAL'
15 include 'COMMON.TORSION'
16 include 'COMMON.FFIELD'
17 include 'COMMON.NAMES'
18 include 'COMMON.SBRIDGE'
19 include 'COMMON.WEIGHTS'
20 include 'COMMON.ENEPS'
21 include 'COMMON.SCCOR'
22 include 'COMMON.SCROT'
25 character*1 onelett(4) /"G","A","P","D"/
27 dimension blower(3,3,maxlob)
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,
36 double precision ip,mp
40 C Set LPRINT=.TRUE. for debugging
41 dwa16=2.0d0**(1.0d0/6.0d0)
44 C Assign virtual-bond length
48 call card_concat(controlcard,.true.)
51 key = wname(i)(:ilen(wname(i)))
52 call reada(controlcard,key(:ilen(key)),ww(i),1.0d0)
54 call reada(controlcard,"D0CM",d0cm,3.78d0)
55 call reada(controlcard,"AKCM",akcm,15.1d0)
56 call reada(controlcard,"AKTH",akth,11.0d0)
57 call reada(controlcard,"AKCT",akct,12.0d0)
58 call reada(controlcard,"V1SS",v1ss,-1.08d0)
59 call reada(controlcard,"V2SS",v2ss,7.61d0)
60 call reada(controlcard,"V3SS",v3ss,13.7d0)
61 call reada(controlcard,"EBR",ebr,-5.50D0)
62 c dyn_ss=(index(controlcard,'DYN_SS').gt.0)
63 write (iout,*) "iparm",iparm," myparm",myparm
65 c dyn_ss_mask(i)=.false.
69 dyn_ssbond_ij(i,j)=1.0d300
72 call reada(controlcard,"HT",Ht,0.0D0)
74 c if(me.eq.king.or..not.out1file) then
75 c print *,'indpdb=',indpdb,' pdbref=',pdbref
77 c If reading not own parameters, skip assignment
78 cc write(iout,*) "KURWA", ww(15)
79 if (iparm.eq.myparm .or. .not.separate_parset) then
82 c Setup weights for UNRES
106 write (iout,*) "wdfa_dist",wdfa_dist," wdfa_tor",wdfa_tor,
107 & " wdfa_nei",wdfa_nei," wdfa_beta",wdfa_beta
110 cc write(iout,*) "KURWA", wstrain,akcm,akth,wsc,dyn_ss
112 call card_concat(controlcard,.false.)
114 c Return if not own parameters
116 if (iparm.ne.myparm .and. separate_parset) return
118 call reads(controlcard,"BONDPAR",bondname_t,bondname)
119 open (ibond,file=bondname_t,status='old')
121 call reads(controlcard,"THETPAR",thetname_t,thetname)
122 open (ithep,file=thetname_t,status='old')
124 call reads(controlcard,"ROTPAR",rotname_t,rotname)
125 open (irotam,file=rotname_t,status='old')
127 call reads(controlcard,"TORPAR",torname_t,torname)
128 open (itorp,file=torname_t,status='old')
130 call reads(controlcard,"TORDPAR",tordname_t,tordname)
131 open (itordp,file=tordname_t,status='old')
133 call reads(controlcard,"SCCORAR",sccorname_t,sccorname)
134 open (isccor,file=sccorname_t,status='old')
136 call reads(controlcard,"FOURIER",fouriername_t,fouriername)
137 open (ifourier,file=fouriername_t,status='old')
139 call reads(controlcard,"ELEPAR",elename_t,elename)
140 open (ielep,file=elename_t,status='old')
142 call reads(controlcard,"SIDEPAR",sidename_t,sidename)
143 open (isidep,file=sidename_t,status='old')
145 call reads(controlcard,"SCPPAR",scpname_t,scpname)
146 open (iscpp,file=scpname_t,status='old')
148 write (iout,*) "Parameter set:",iparm
149 write (iout,*) "Energy-term weights:"
151 write (iout,'(a16,f10.5)') wname(i),ww(i)
153 write (iout,*) "Sidechain potential file : ",
154 & sidename_t(:ilen(sidename_t))
156 write (iout,*) "SCp potential file : ",
157 & scpname_t(:ilen(scpname_t))
159 write (iout,*) "Electrostatic potential file : ",
160 & elename_t(:ilen(elename_t))
161 write (iout,*) "Cumulant coefficient file : ",
162 & fouriername_t(:ilen(fouriername_t))
163 write (iout,*) "Torsional parameter file : ",
164 & torname_t(:ilen(torname_t))
165 write (iout,*) "Double torsional parameter file : ",
166 & tordname_t(:ilen(tordname_t))
167 write (iout,*) "Backbone-rotamer parameter file : ",
168 & sccorname(:ilen(sccorname))
169 write (iout,*) "Bond & inertia constant file : ",
170 & bondname_t(:ilen(bondname_t))
171 write (iout,*) "Bending parameter file : ",
172 & thetname_t(:ilen(thetname_t))
173 write (iout,*) "Rotamer parameter file : ",
174 & rotname_t(:ilen(rotname_t))
177 c Read the virtual-bond parameters, masses, and moments of inertia
178 c and Stokes' radii of the peptide group and side chains
181 read (ibond,*) vbldp0,akp
184 read (ibond,*) vbldsc0(1,i),aksc(1,i)
185 dsc(i) = vbldsc0(1,i)
189 dsc_inv(i)=1.0D0/dsc(i)
193 read (ibond,*) ijunk,vbldp0,akp,rjunk
195 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
197 dsc(i) = vbldsc0(1,i)
201 dsc_inv(i)=1.0D0/dsc(i)
206 write(iout,'(/a/)')"Force constants virtual bonds:"
207 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
209 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
211 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
212 & vbldsc0(1,i),aksc(1,i),abond0(1,i)
214 write (iout,'(13x,3f10.5)')
215 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
221 C Read the parameters of the probability distribution/energy expression
222 C of the virtual-bond valence angles theta
225 read (ithep,*) a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
226 read (ithep,*) (polthet(j,i),j=0,3)
227 read (ithep,*) (gthet(j,i),j=1,3)
228 read (ithep,*) theta0(i),sig0(i),sigc0(i)
234 c & 'Parameters of the virtual-bond valence angles:'
235 c write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
236 c & ' ATHETA0 ',' A1 ',' A2 ',
239 c write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
240 c & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
242 c write (iout,'(/a/9x,5a/79(1h-))')
243 c & 'Parameters of the expression for sigma(theta_c):',
244 c & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
245 c & ' ALPH3 ',' SIGMA0C '
247 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
248 c & (polthet(j,i),j=0,3),sigc0(i)
250 c write (iout,'(/a/9x,5a/79(1h-))')
251 c & 'Parameters of the second gaussian:',
252 c & ' THETA0 ',' SIGMA0 ',' G1 ',
255 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
256 c & sig0(i),(gthet(j,i),j=1,3)
259 & 'Parameters of the virtual-bond valence angles:'
260 write (iout,'(/a/9x,5a/79(1h-))')
261 & 'Coefficients of expansion',
262 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
263 & ' b1*10^1 ',' b2*10^1 '
265 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
266 & a0thet(i),(100*athet(j,i),j=1,2),(10*bthet(j,i),j=1,2)
268 write (iout,'(/a/9x,5a/79(1h-))')
269 & 'Parameters of the expression for sigma(theta_c):',
270 & ' alpha0 ',' alph1 ',' alph2 ',
271 & ' alhp3 ',' sigma0c '
273 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
274 & (polthet(j,i),j=0,3),sigc0(i)
276 write (iout,'(/a/9x,5a/79(1h-))')
277 & 'Parameters of the second gaussian:',
278 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
281 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
282 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
287 C Read the parameters of Utheta determined from ab initio surfaces
288 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
290 read (ithep,*) nthetyp,ntheterm,ntheterm2,
291 & ntheterm3,nsingle,ndouble
292 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
293 read (ithep,*) (ithetyp(i),i=1,ntyp1)
299 aathet(l,i,j,k)=0.0d0
303 bbthet(m,l,i,j,k)=0.0d0
304 ccthet(m,l,i,j,k)=0.0d0
305 ddthet(m,l,i,j,k)=0.0d0
306 eethet(m,l,i,j,k)=0.0d0
312 ffthet(mm,m,l,i,j,k)=0.0d0
313 ggthet(mm,m,l,i,j,k)=0.0d0
323 read (ithep,'(3a)') res1,res2,res3
324 read (ithep,*) aa0thet(i,j,k)
325 read (ithep,*)(aathet(l,i,j,k),l=1,ntheterm)
327 & ((bbthet(lll,ll,i,j,k),lll=1,nsingle),
328 & (ccthet(lll,ll,i,j,k),lll=1,nsingle),
329 & (ddthet(lll,ll,i,j,k),lll=1,nsingle),
330 & (eethet(lll,ll,i,j,k),lll=1,nsingle),ll=1,ntheterm2)
332 & (((ffthet(llll,lll,ll,i,j,k),ffthet(lll,llll,ll,i,j,k),
333 & ggthet(llll,lll,ll,i,j,k),ggthet(lll,llll,ll,i,j,k),
334 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
339 C For dummy ends assign glycine-type coefficients of theta-only terms; the
340 C coefficients of theta-and-gamma-dependent terms are zero.
345 aathet(l,i,j,nthetyp+1)=aathet(l,i,j,1)
346 aathet(l,nthetyp+1,i,j)=aathet(l,1,i,j)
348 aa0thet(i,j,nthetyp+1)=aa0thet(i,j,1)
349 aa0thet(nthetyp+1,i,j)=aa0thet(1,i,j)
352 aathet(l,nthetyp+1,i,nthetyp+1)=aathet(l,1,i,1)
354 aa0thet(nthetyp+1,i,nthetyp+1)=aa0thet(1,i,1)
357 C Control printout of the coefficients of virtual-bond-angle potentials
360 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
364 write (iout,'(//4a)')
365 & 'Type ',onelett(i),onelett(j),onelett(k)
366 write (iout,'(//a,10x,a)') " l","a[l]"
367 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k)
368 write (iout,'(i2,1pe15.5)')
369 & (l,aathet(l,i,j,k),l=1,ntheterm)
371 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
372 & "b",l,"c",l,"d",l,"e",l
374 write (iout,'(i2,4(1pe15.5))') m,
375 & bbthet(m,l,i,j,k),ccthet(m,l,i,j,k),
376 & ddthet(m,l,i,j,k),eethet(m,l,i,j,k)
380 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
381 & "f+",l,"f-",l,"g+",l,"g-",l
384 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
385 & ffthet(n,m,l,i,j,k),ffthet(m,n,l,i,j,k),
386 & ggthet(n,m,l,i,j,k),ggthet(m,n,l,i,j,k)
399 C Read the parameters of the probability distribution/energy expression
400 C of the side chains.
403 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
407 dsc_inv(i)=1.0D0/dsc(i)
418 read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
420 read (irotam,*) bsc(j,i)
421 read (irotam,*) (censc(k,j,i),k=1,3),
422 & ((blower(k,l,j),l=1,k),k=1,3)
429 akl=akl+blower(k,m,j)*blower(l,m,j)
440 write (iout,'(/a)') 'Parameters of side-chain local geometry'
444 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
445 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
446 c write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
447 c write (iout,'(a,f10.4,4(16x,f10.4))')
448 c & 'Center ',(bsc(j,i),j=1,nlobi)
449 c write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
450 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
451 & 'log h',(bsc(j,i),j=1,nlobi)
452 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
453 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
460 c blower(k,l,j)=gaussc(ind,j,i)
465 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
466 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
473 C Read scrot parameters for potentials determined from all-atom AM1 calculations
474 C added by Urszula Kozlowska 07/11/2007
482 read(irotam,*) sc_parmin(j,i)
490 C Read torsional parameters in old format
492 read (itorp,*) ntortyp,nterm_old
493 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
494 read (itorp,*) (itortyp(i),i=1,ntyp)
499 read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
505 write (iout,'(/a/)') 'Torsional constants:'
508 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
509 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
517 C Read torsional parameters
519 read (itorp,*) ntortyp
520 read (itorp,*) (itortyp(i),i=1,ntyp)
521 write (iout,*) 'ntortyp',ntortyp
524 read (itorp,*) nterm(i,j),nlor(i,j)
528 read (itorp,*) kk,v1(k,i,j),v2(k,i,j)
529 v0ij=v0ij+si*v1(k,i,j)
533 read (itorp,*) kk,vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
534 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
541 write (iout,'(/a/)') 'Torsional constants:'
544 write (iout,*) 'ityp',i,' jtyp',j
545 write (iout,*) 'Fourier constants'
547 write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
549 write (iout,*) 'Lorenz constants'
551 write (iout,'(3(1pe15.5))')
552 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
558 C 6/23/01 Read parameters for double torsionals
563 read (itordp,'(3a1)') t1,t2,t3
564 if (t1.ne.onelett(i) .or. t2.ne.onelett(j)
565 & .or. t3.ne.onelett(k)) then
566 write (iout,*) "Error in double torsional parameter file",
568 stop "Error in double torsional parameter file"
570 read (itordp,*) ntermd_1(i,j,k),ntermd_2(i,j,k)
571 read (itordp,*) (v1c(1,l,i,j,k),l=1,ntermd_1(i,j,k))
572 read (itordp,*) (v1s(1,l,i,j,k),l=1,ntermd_1(i,j,k))
573 read (itordp,*) (v1c(2,l,i,j,k),l=1,ntermd_1(i,j,k))
574 read (itordp,*) (v1s(2,l,i,j,k),l=1,ntermd_1(i,j,k))
575 read (itordp,*) ((v2c(l,m,i,j,k),v2c(m,l,i,j,k),
576 & v2s(l,m,i,j,k),v2s(m,l,i,j,k),m=1,l-1),l=1,ntermd_2(i,j,k))
582 write (iout,*) 'Constants for double torsionals'
586 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
587 & ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k)
589 write (iout,*) 'Single angles:'
590 do l=1,ntermd_1(i,j,k)
591 write (iout,'(i5,2f10.5,5x,2f10.5)') l,
592 & v1c(1,l,i,j,k),v1s(1,l,i,j,k),
593 & v1c(2,l,i,j,k),v1s(2,l,i,j,k)
596 write (iout,*) 'Pairs of angles:'
597 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
598 do l=1,ntermd_2(i,j,k)
599 write (iout,'(i5,20f10.5)')
600 & l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k))
603 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
604 do l=1,ntermd_2(i,j,k)
605 write (iout,'(i5,20f10.5)')
606 & l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k))
614 C Read of Side-chain backbone correlation parameters
615 C Modified 11 May 2012 by Adasko
618 read (isccor,*) nsccortyp
619 read (isccor,*) (isccortyp(i),i=1,ntyp)
620 c write (iout,*) 'ntortyp',ntortyp
622 cc maxinter is maximum interaction sites
626 read (isccor,*) nterm_sccor(i,j),nlor_sccor(i,j)
630 do k=1,nterm_sccor(i,j)
631 read (isccor,*) kk,v1sccor(k,l,i,j)
633 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
636 do k=1,nlor_sccor(i,j)
637 read (isccor,*) kk,vlor1sccor(k,i,j),
638 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
639 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
640 &(1+vlor3sccor(k,i,j)**2)
642 v0sccor(i,j)=v0ijsccor
649 write (iout,'(/a/)') 'Torsional constants:'
652 write (iout,*) 'ityp',i,' jtyp',j
653 write (iout,*) 'Fourier constants'
654 do k=1,nterm_sccor(i,j)
655 write (iout,'(2(1pe15.5))')v1sccor(k,l,i,j),v2sccor(k,l,i,j)
657 write (iout,*) 'Lorenz constants'
658 do k=1,nlor_sccor(i,j)
659 write (iout,'(3(1pe15.5))')
660 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
667 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
668 C interaction energy of the Gly, Ala, and Pro prototypes.
670 read (ifourier,*) nloctyp
673 read (ifourier,*) (b(ii,i),ii=1,13)
675 write (iout,*) 'Type',i
676 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
680 B1tilde(1,i) = b(3,i)
681 B1tilde(2,i) =-b(5,i)
690 Ctilde(2,1,i)=-b(9,i)
698 Dtilde(2,1,i)=-b(8,i)
700 EE(1,1,i)= b(10,i)+b(11,i)
701 EE(2,2,i)=-b(10,i)+b(11,i)
702 EE(2,1,i)= b(12,i)-b(13,i)
703 EE(1,2,i)= b(12,i)+b(13,i)
707 write (iout,*) 'Type',i
709 c write (iout,'(f10.5)') B1(:,i)
710 write(iout,*) B1(1,i),B1(2,i)
712 c write (iout,'(f10.5)') B2(:,i)
713 write(iout,*) B2(1,i),B2(2,i)
716 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
720 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
724 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
729 C Read electrostatic-interaction parameters
732 write (iout,'(/a)') 'Electrostatic interaction constants:'
733 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
734 & 'IT','JT','APP','BPP','AEL6','AEL3'
736 read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
737 read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
738 read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
739 read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
744 app (i,j)=epp(i,j)*rri*rri
745 bpp (i,j)=-2.0D0*epp(i,j)*rri
746 ael6(i,j)=elpp6(i,j)*4.2D0**6
747 ael3(i,j)=elpp3(i,j)*4.2D0**3
748 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
749 & ael6(i,j),ael3(i,j)
753 C Read side-chain interaction parameters.
755 read (isidep,*) ipot,expon
756 if (ipot.lt.1 .or. ipot.gt.5) then
757 write (iout,'(2a)') 'Error while reading SC interaction',
758 & 'potential file - unknown potential type.'
762 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
763 & ', exponents are ',expon,2*expon
764 goto (10,20,30,30,40) ipot
765 C----------------------- LJ potential ---------------------------------
766 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
768 write (iout,'(/a/)') 'Parameters of the LJ potential:'
769 write (iout,'(a/)') 'The epsilon array:'
770 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
771 write (iout,'(/a)') 'One-body parameters:'
772 write (iout,'(a,4x,a)') 'residue','sigma'
773 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
776 C----------------------- LJK potential --------------------------------
777 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
778 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
780 write (iout,'(/a/)') 'Parameters of the LJK potential:'
781 write (iout,'(a/)') 'The epsilon array:'
782 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
783 write (iout,'(/a)') 'One-body parameters:'
784 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
785 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
789 C---------------------- GB or BP potential -----------------------------
790 30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
791 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
793 C For the GB potential convert sigma'**2 into chi'
796 chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
800 write (iout,'(/a/)') 'Parameters of the BP potential:'
801 write (iout,'(a/)') 'The epsilon array:'
802 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
803 write (iout,'(/a)') 'One-body parameters:'
804 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
806 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
807 & chip(i),alp(i),i=1,ntyp)
810 C--------------------- GBV potential -----------------------------------
811 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
812 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
813 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
815 write (iout,'(/a/)') 'Parameters of the GBV potential:'
816 write (iout,'(a/)') 'The epsilon array:'
817 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
818 write (iout,'(/a)') 'One-body parameters:'
819 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
820 & 's||/s_|_^2',' chip ',' alph '
821 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
822 & sigii(i),chip(i),alp(i),i=1,ntyp)
826 C-----------------------------------------------------------------------
827 C Calculate the "working" parameters of SC interactions.
835 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
836 sigma(j,i)=sigma(i,j)
837 rs0(i,j)=dwa16*sigma(i,j)
841 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
842 & 'Working parameters of the SC interactions:',
843 & ' a ',' b ',' augm ',' sigma ',' r0 ',
848 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
857 sigeps=dsign(1.0D0,epsij)
859 aa(i,j)=epsij*rrij*rrij
860 bb(i,j)=-sigeps*epsij*rrij
868 ratsig1=sigt2sq/sigt1sq
869 ratsig2=1.0D0/ratsig1
870 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
871 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
872 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
876 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
877 sigmaii(i,j)=rsum_max
878 sigmaii(j,i)=rsum_max
880 c sigmaii(i,j)=r0(i,j)
881 c sigmaii(j,i)=r0(i,j)
883 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
884 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
885 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
886 augm(i,j)=epsij*r_augm**(2*expon)
887 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
894 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
895 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
896 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
901 C Define the SC-p interaction constants and SS bond potentials
904 ss_depth=ebr/wsc-0.25*eps(1,1)
905 Ht=Ht/wsc-0.25*eps(1,1)
906 akcm=akcm*wstrain/wsc
907 akth=akth*wstrain/wsc
908 akct=akct*wstrain/wsc
909 v1ss=v1ss*wstrain/wsc
910 v2ss=v2ss*wstrain/wsc
911 v3ss=v3ss*wstrain/wsc
913 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
915 write (iout,*) "Parameters of the SS-bond potential:"
916 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
918 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
919 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
920 write (iout,*)" HT",Ht
924 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
926 c aad(i,1)=0.3D0*4.0D0**12
927 C Following line for constants currently implemented
928 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
929 aad(i,1)=1.5D0*4.0D0**12
930 c aad(i,1)=0.17D0*5.6D0**12
932 C "Soft" SC-p repulsion
934 C Following line for constants currently implemented
935 c aad(i,1)=0.3D0*4.0D0**6
936 C "Hard" SC-p repulsion
937 bad(i,1)=3.0D0*4.0D0**6
938 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
947 C 8/9/01 Read the SC-p interaction constants from file
950 read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
953 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
954 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
955 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
956 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
960 write (iout,*) "Parameters of SC-p interactions:"
962 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
963 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
968 C Define the constants of the disulfide bridge
972 c Old arbitrary potential - commented out.
977 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
978 c energy surface of diethyl disulfide.
979 c A. Liwo and U. Kozlowska, 11/24/03
990 c write (iout,'(/a)') "Disulfide bridge parameters:"
991 c write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
992 c write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
993 c write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
994 c write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,