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)
55 write (iout,*) "iparm",iparm," myparm",myparm
56 c If reading not own parameters, skip assignment
58 if (iparm.eq.myparm .or. .not.separate_parset) then
61 c Setup weights for UNRES
84 call card_concat(controlcard,.false.)
86 c Return if not own parameters
88 if (iparm.ne.myparm .and. separate_parset) return
90 call reads(controlcard,"BONDPAR",bondname_t,bondname)
91 open (ibond,file=bondname_t,status='old')
93 call reads(controlcard,"THETPAR",thetname_t,thetname)
94 open (ithep,file=thetname_t,status='old')
96 call reads(controlcard,"ROTPAR",rotname_t,rotname)
97 open (irotam,file=rotname_t,status='old')
99 call reads(controlcard,"TORPAR",torname_t,torname)
100 open (itorp,file=torname_t,status='old')
102 call reads(controlcard,"TORDPAR",tordname_t,tordname)
103 open (itordp,file=tordname_t,status='old')
105 call reads(controlcard,"SCCORAR",sccorname_t,sccorname)
106 open (isccor,file=sccorname_t,status='old')
108 call reads(controlcard,"FOURIER",fouriername_t,fouriername)
109 open (ifourier,file=fouriername_t,status='old')
111 call reads(controlcard,"ELEPAR",elename_t,elename)
112 open (ielep,file=elename_t,status='old')
114 call reads(controlcard,"SIDEPAR",sidename_t,sidename)
115 open (isidep,file=sidename_t,status='old')
117 call reads(controlcard,"SCPPAR",scpname_t,scpname)
118 open (iscpp,file=scpname_t,status='old')
120 write (iout,*) "Parameter set:",iparm
121 write (iout,*) "Energy-term weights:"
123 write (iout,'(a16,f10.5)') wname(i),ww(i)
125 write (iout,*) "Sidechain potential file : ",
126 & sidename_t(:ilen(sidename_t))
128 write (iout,*) "SCp potential file : ",
129 & scpname_t(:ilen(scpname_t))
131 write (iout,*) "Electrostatic potential file : ",
132 & elename_t(:ilen(elename_t))
133 write (iout,*) "Cumulant coefficient file : ",
134 & fouriername_t(:ilen(fouriername_t))
135 write (iout,*) "Torsional parameter file : ",
136 & torname_t(:ilen(torname_t))
137 write (iout,*) "Double torsional parameter file : ",
138 & tordname_t(:ilen(tordname_t))
139 write (iout,*) "Backbone-rotamer parameter file : ",
140 & sccorname(:ilen(sccorname))
141 write (iout,*) "Bond & inertia constant file : ",
142 & bondname_t(:ilen(bondname_t))
143 write (iout,*) "Bending parameter file : ",
144 & thetname_t(:ilen(thetname_t))
145 write (iout,*) "Rotamer parameter file : ",
146 & rotname_t(:ilen(rotname_t))
149 c Read the virtual-bond parameters, masses, and moments of inertia
150 c and Stokes' radii of the peptide group and side chains
153 read (ibond,*,end=110,err=110) vbldp0,akp
156 read (ibond,*,end=110,err=110) vbldsc0(1,i),aksc(1,i)
157 dsc(i) = vbldsc0(1,i)
161 dsc_inv(i)=1.0D0/dsc(i)
165 read (ibond,*,end=110,err=110) ijunk,vbldp0,akp,rjunk
167 read (ibond,*,end=110,err=110) nbondterm(i),(vbldsc0(j,i),
168 & aksc(j,i),abond0(j,i),
170 dsc(i) = vbldsc0(1,i)
174 dsc_inv(i)=1.0D0/dsc(i)
179 write(iout,'(/a/)')"Force constants virtual bonds:"
180 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
182 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
184 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
185 & vbldsc0(1,i),aksc(1,i),abond0(1,i)
187 write (iout,'(13x,3f10.5)')
188 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
194 C Read the parameters of the probability distribution/energy expression
195 C of the virtual-bond valence angles theta
198 read (ithep,*,end=111,err=111) a0thet(i),(athet(j,i),j=1,2),
200 read (ithep,*,end=111,err=111) (polthet(j,i),j=0,3)
201 read (ithep,*,end=111,err=111) (gthet(j,i),j=1,3)
202 read (ithep,*,end=111,err=111) theta0(i),sig0(i),sigc0(i)
208 c & 'Parameters of the virtual-bond valence angles:'
209 c write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
210 c & ' ATHETA0 ',' A1 ',' A2 ',
213 c write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
214 c & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
216 c write (iout,'(/a/9x,5a/79(1h-))')
217 c & 'Parameters of the expression for sigma(theta_c):',
218 c & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
219 c & ' ALPH3 ',' SIGMA0C '
221 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
222 c & (polthet(j,i),j=0,3),sigc0(i)
224 c write (iout,'(/a/9x,5a/79(1h-))')
225 c & 'Parameters of the second gaussian:',
226 c & ' THETA0 ',' SIGMA0 ',' G1 ',
229 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
230 c & sig0(i),(gthet(j,i),j=1,3)
233 & 'Parameters of the virtual-bond valence angles:'
234 write (iout,'(/a/9x,5a/79(1h-))')
235 & 'Coefficients of expansion',
236 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
237 & ' b1*10^1 ',' b2*10^1 '
239 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
240 & a0thet(i),(100*athet(j,i),j=1,2),(10*bthet(j,i),j=1,2)
242 write (iout,'(/a/9x,5a/79(1h-))')
243 & 'Parameters of the expression for sigma(theta_c):',
244 & ' alpha0 ',' alph1 ',' alph2 ',
245 & ' alhp3 ',' sigma0c '
247 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
248 & (polthet(j,i),j=0,3),sigc0(i)
250 write (iout,'(/a/9x,5a/79(1h-))')
251 & 'Parameters of the second gaussian:',
252 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
255 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
256 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
261 C Read the parameters of Utheta determined from ab initio surfaces
262 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
264 read (ithep,*,end=111,err=111) nthetyp,ntheterm,ntheterm2,
265 & ntheterm3,nsingle,ndouble
266 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
267 read (ithep,*,end=111,err=111) (ithetyp(i),i=1,ntyp1)
273 aathet(l,i,j,k)=0.0d0
277 bbthet(m,l,i,j,k)=0.0d0
278 ccthet(m,l,i,j,k)=0.0d0
279 ddthet(m,l,i,j,k)=0.0d0
280 eethet(m,l,i,j,k)=0.0d0
286 ffthet(mm,m,l,i,j,k)=0.0d0
287 ggthet(mm,m,l,i,j,k)=0.0d0
297 read (ithep,'(3a)',end=111,err=111) res1,res2,res3
298 read (ithep,*,end=111,err=111) aa0thet(i,j,k)
299 read (ithep,*,end=111,err=111)(aathet(l,i,j,k),l=1,ntheterm)
300 read (ithep,*,end=111,err=111)
301 & ((bbthet(lll,ll,i,j,k),lll=1,nsingle),
302 & (ccthet(lll,ll,i,j,k),lll=1,nsingle),
303 & (ddthet(lll,ll,i,j,k),lll=1,nsingle),
304 & (eethet(lll,ll,i,j,k),lll=1,nsingle),ll=1,ntheterm2)
305 read (ithep,*,end=111,err=111)
306 & (((ffthet(llll,lll,ll,i,j,k),ffthet(lll,llll,ll,i,j,k),
307 & ggthet(llll,lll,ll,i,j,k),ggthet(lll,llll,ll,i,j,k),
308 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
313 C For dummy ends assign glycine-type coefficients of theta-only terms; the
314 C coefficients of theta-and-gamma-dependent terms are zero.
319 aathet(l,i,j,nthetyp+1)=aathet(l,i,j,1)
320 aathet(l,nthetyp+1,i,j)=aathet(l,1,i,j)
322 aa0thet(i,j,nthetyp+1)=aa0thet(i,j,1)
323 aa0thet(nthetyp+1,i,j)=aa0thet(1,i,j)
326 aathet(l,nthetyp+1,i,nthetyp+1)=aathet(l,1,i,1)
328 aa0thet(nthetyp+1,i,nthetyp+1)=aa0thet(1,i,1)
331 C Control printout of the coefficients of virtual-bond-angle potentials
334 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
338 write (iout,'(//4a)')
339 & 'Type ',onelett(i),onelett(j),onelett(k)
340 write (iout,'(//a,10x,a)') " l","a[l]"
341 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k)
342 write (iout,'(i2,1pe15.5)')
343 & (l,aathet(l,i,j,k),l=1,ntheterm)
345 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
346 & "b",l,"c",l,"d",l,"e",l
348 write (iout,'(i2,4(1pe15.5))') m,
349 & bbthet(m,l,i,j,k),ccthet(m,l,i,j,k),
350 & ddthet(m,l,i,j,k),eethet(m,l,i,j,k)
354 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
355 & "f+",l,"f-",l,"g+",l,"g-",l
358 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
359 & ffthet(n,m,l,i,j,k),ffthet(m,n,l,i,j,k),
360 & ggthet(n,m,l,i,j,k),ggthet(m,n,l,i,j,k)
373 C Read the parameters of the probability distribution/energy expression
374 C of the side chains.
377 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
381 dsc_inv(i)=1.0D0/dsc(i)
392 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
393 & ((blower(k,l,1),l=1,k),k=1,3)
395 read (irotam,*,end=112,err=112) bsc(j,i)
396 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
397 & ((blower(k,l,j),l=1,k),k=1,3)
404 akl=akl+blower(k,m,j)*blower(l,m,j)
415 write (iout,'(/a)') 'Parameters of side-chain local geometry'
419 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
420 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
421 c write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
422 c write (iout,'(a,f10.4,4(16x,f10.4))')
423 c & 'Center ',(bsc(j,i),j=1,nlobi)
424 c write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
425 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
426 & 'log h',(bsc(j,i),j=1,nlobi)
427 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
428 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
435 c blower(k,l,j)=gaussc(ind,j,i)
440 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
441 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
448 C Read scrot parameters for potentials determined from all-atom AM1 calculations
449 C added by Urszula Kozlowska 07/11/2007
452 read (irotam,*,end=112,err=112)
454 read (irotam,*,end=112,err=112)
457 read(irotam,*,end=112,err=112) sc_parmin(j,i)
465 C Read torsional parameters in old format
467 read (itorp,*,end=113,err=113) ntortyp,nterm_old
468 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
469 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
472 read (itorp,'(a)',end=113,err=113)
474 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
480 write (iout,'(/a/)') 'Torsional constants:'
483 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
484 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
492 C Read torsional parameters
494 read (itorp,*,end=113,err=113) ntortyp
495 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
496 write (iout,*) 'ntortyp',ntortyp
499 read (itorp,*,end=113,err=113) nterm(i,j),nlor(i,j)
503 read (itorp,*,end=113,err=113) kk,v1(k,i,j),v2(k,i,j)
504 v0ij=v0ij+si*v1(k,i,j)
508 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),vlor2(k,i,j),
510 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
517 write (iout,'(/a/)') 'Torsional constants:'
520 write (iout,*) 'ityp',i,' jtyp',j
521 write (iout,*) 'Fourier constants'
523 write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
525 write (iout,*) 'Lorenz constants'
527 write (iout,'(3(1pe15.5))')
528 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
534 C 6/23/01 Read parameters for double torsionals
539 read (itordp,'(3a1)',end=112,err=112) t1,t2,t3
540 if (t1.ne.onelett(i) .or. t2.ne.onelett(j)
541 & .or. t3.ne.onelett(k)) then
542 write (iout,*) "Error in double torsional parameter file",
544 stop "Error in double torsional parameter file"
546 read (itordp,*,end=114,err=114) ntermd_1(i,j,k),
548 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k),
549 & l=1,ntermd_1(i,j,k))
550 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k),l=1,
552 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k),
553 & l=1,ntermd_1(i,j,k))
554 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k),l=1,
556 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k),
558 & v2s(l,m,i,j,k),v2s(m,l,i,j,k),m=1,l-1),l=1,ntermd_2(i,j,k))
564 write (iout,*) 'Constants for double torsionals'
568 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
569 & ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k)
571 write (iout,*) 'Single angles:'
572 do l=1,ntermd_1(i,j,k)
573 write (iout,'(i5,2f10.5,5x,2f10.5)') l,
574 & v1c(1,l,i,j,k),v1s(1,l,i,j,k),
575 & v1c(2,l,i,j,k),v1s(2,l,i,j,k)
578 write (iout,*) 'Pairs of angles:'
579 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
580 do l=1,ntermd_2(i,j,k)
581 write (iout,'(i5,20f10.5)')
582 & l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k))
585 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
586 do l=1,ntermd_2(i,j,k)
587 write (iout,'(i5,20f10.5)')
588 & l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k))
596 C Read of Side-chain backbone correlation parameters
597 C Modified 11 May 2012 by Adasko
600 read (isccor,*,end=119,err=119) nsccortyp
601 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
602 c write (iout,*) 'ntortyp',ntortyp
604 cc maxinter is maximum interaction sites
608 read (isccor,*,end=119,err=119) nterm_sccor(i,j),
612 do k=1,nterm_sccor(i,j)
613 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
615 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
618 do k=1,nlor_sccor(i,j)
619 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
620 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
621 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
622 &(1+vlor3sccor(k,i,j)**2)
624 v0sccor(i,j)=v0ijsccor
631 write (iout,'(/a/)') 'Torsional constants:'
634 write (iout,*) 'ityp',i,' jtyp',j
635 write (iout,*) 'Fourier constants'
636 do k=1,nterm_sccor(i,j)
637 write(iout,'(2(1pe15.5))')(v1sccor(k,l,i,j),v2sccor(k,l,i,j),
640 write (iout,*) 'Lorenz constants'
641 do k=1,nlor_sccor(i,j)
642 write (iout,'(3(1pe15.5))')
643 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
650 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
651 C interaction energy of the Gly, Ala, and Pro prototypes.
653 read (ifourier,*,end=115,err=115) nloctyp
655 read (ifourier,*,end=115,err=115)
656 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
658 read (ifourier,*,end=115,err=115) (bnew1(ii,1,i),ii=1,3)
659 read (ifourier,*,end=115,err=115) (bnew2(ii,1,i),ii=1,3)
660 read (ifourier,*,end=115,err=115) (bnew1(ii,2,i),ii=1,1)
661 read (ifourier,*,end=115,err=115) (bnew2(ii,2,i),ii=1,1)
662 read (ifourier,*,end=115,err=115) (eenew(ii,i),ii=1,1)
665 write (iout,*) 'Type',i
666 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
671 B1tilde(1,i) = b(3,i)
672 B1tilde(2,i) =-b(5,i)
682 Ctilde(2,1,i)=-b(9,i)
690 Dtilde(2,1,i)=-b(8,i)
693 EEold(1,1,i)= b(10,i)+b(11,i)
694 EEold(2,2,i)=-b(10,i)+b(11,i)
695 EEold(2,1,i)= b(12,i)-b(13,i)
696 EEold(1,2,i)= b(12,i)+b(13,i)
697 EEold(1,1,-i)= b(10,i)+b(11,i)
698 EEold(2,2,-i)=-b(10,i)+b(11,i)
699 EEold(2,1,-i)=-b(12,i)+b(13,i)
700 EEold(1,2,-i)=-b(12,i)-b(13,i)
702 EE(1,1,i)= b(10,i)+b(11,i)
703 EE(2,2,i)=-b(10,i)+b(11,i)
704 EE(2,1,i)= b(12,i)-b(13,i)
705 EE(1,2,i)= b(12,i)+b(13,i)
710 write (iout,*) 'Type',i
712 c write (iout,'(f10.5)') B1(:,i)
713 write(iout,*) B1(1,i),B1(2,i)
715 c write (iout,'(f10.5)') B2(:,i)
716 write(iout,*) B2(1,i),B2(2,i)
719 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
723 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
727 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
732 C Read electrostatic-interaction parameters
735 write (iout,'(/a)') 'Electrostatic interaction constants:'
736 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
737 & 'IT','JT','APP','BPP','AEL6','AEL3'
739 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
740 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
741 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
742 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
747 app (i,j)=epp(i,j)*rri*rri
748 bpp (i,j)=-2.0D0*epp(i,j)*rri
749 ael6(i,j)=elpp6(i,j)*4.2D0**6
750 ael3(i,j)=elpp3(i,j)*4.2D0**3
751 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
752 & ael6(i,j),ael3(i,j)
756 C Read side-chain interaction parameters.
758 read (isidep,*,end=117,err=117) ipot,expon
759 if (ipot.lt.1 .or. ipot.gt.6) then
760 write (iout,'(2a)') 'Error while reading SC interaction',
761 & 'potential file - unknown potential type.'
765 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
766 & ', exponents are ',expon,2*expon
767 goto (10,20,30,30,40,50) ipot
768 C----------------------- LJ potential ---------------------------------
769 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
770 & (sigma0(i),i=1,ntyp)
772 write (iout,'(/a/)') 'Parameters of the LJ potential:'
773 write (iout,'(a/)') 'The epsilon array:'
774 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
775 write (iout,'(/a)') 'One-body parameters:'
776 write (iout,'(a,4x,a)') 'residue','sigma'
777 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
780 C----------------------- LJK potential --------------------------------
781 20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
782 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
784 write (iout,'(/a/)') 'Parameters of the LJK potential:'
785 write (iout,'(a/)') 'The epsilon array:'
786 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
787 write (iout,'(/a)') 'One-body parameters:'
788 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
789 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
793 C---------------------- GB or BP potential -----------------------------
794 30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
795 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
797 C For the GB potential convert sigma'**2 into chi'
800 chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
804 write (iout,'(/a/)') 'Parameters of the BP potential:'
805 write (iout,'(a/)') 'The epsilon array:'
806 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
807 write (iout,'(/a)') 'One-body parameters:'
808 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
810 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
811 & chip(i),alp(i),i=1,ntyp)
814 C--------------------- GBV potential -----------------------------------
815 40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
816 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
817 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
819 write (iout,'(/a/)') 'Parameters of the GBV potential:'
820 write (iout,'(a/)') 'The epsilon array:'
821 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
822 write (iout,'(/a)') 'One-body parameters:'
823 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
824 & 's||/s_|_^2',' chip ',' alph '
825 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
826 & sigii(i),chip(i),alp(i),i=1,ntyp)
829 C--------------------- Momo potential -----------------------------------
833 read (isidep,*,end=116,err=116) (icharge(i),i=1,ntyp)
834 c write (2,*) "icharge",(icharge(i),i=1,ntyp)
837 c! write (*,*) "Im in ", i, " ", j
838 read(isidep,*,end=116,err=116)
839 & eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i),
840 & (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j),
841 & chis(i,j),chis(j,i),
842 & nstate(i,j),(wstate(k,i,j),k=1,4),
847 & dtail(1,i,j),dtail(2,i,j),
848 & epshead(i,j),sig0head(i,j),
849 & rborn(i,j),rborn(j,i),
850 & (wqdip(k,i,j),k=1,2),wquad(i,j),
851 & alphapol(i,j),alphapol(j,i),
852 & (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j)
855 c! write (*,*) "nstate gly-gly", nstate(10,10)
856 c! THIS LOOP FILLS PARAMETERS FOR PAIRS OF AA's NOT EXPLICITLY
857 c! DEFINED IN SCPARM.MOMO. IT DOES SO BY TAKING THEM FROM SYMMETRIC
858 c! PAIR, FOR EG. IF ARG-HIS IS BLANK, IT WILL TAKE PARAMETERS
861 c! DISABLE IT AT >>YOUR OWN PERIL<<
866 sigma(i,j) = sigma(j,i)
867 nstate(i,j) = nstate(j,i)
868 sigmap1(i,j) = sigmap1(j,i)
869 sigmap2(i,j) = sigmap2(j,i)
870 sigiso1(i,j) = sigiso1(j,i)
871 sigiso2(i,j) = sigiso2(j,i)
874 alphasur(k,i,j) = alphasur(k,j,i)
875 wstate(k,i,j) = wstate(k,j,i)
876 alphiso(k,i,j) = alphiso(k,j,i)
879 dhead(2,1,i,j) = dhead(1,1,j,i)
880 dhead(2,2,i,j) = dhead(1,2,j,i)
881 dhead(1,1,i,j) = dhead(2,1,j,i)
882 dhead(1,2,i,j) = dhead(2,2,j,i)
883 dtail(1,i,j) = dtail(1,j,i)
884 dtail(2,i,j) = dtail(2,j,i)
887 c! write (*,*) "dhead(k,l,j,i) = ", dhead(k,l,j,i)
888 c! write (*,*) "dhead(k,l,i,j) = ", dhead(k,l,i,j)
889 c! dhead(l,k,i,j) = dhead(k,l,j,i)
893 epshead(i,j) = epshead(j,i)
894 sig0head(i,j) = sig0head(j,i)
897 wqdip(k,i,j) = wqdip(k,j,i)
900 wquad(i,j) = wquad(j,i)
901 epsintab(i,j) = epsintab(j,i)
906 if (.not.lprint) goto 70
908 & "Parameters of the new physics-based SC-SC interaction potential"
909 write (iout,'(/7a)') 'Residues',' epsGB',' rGB',
910 & ' chi1GB',' chi2GB',' chip1GB',' chip2GB'
913 write (iout,'(2(a3,1x),1pe10.3,5(0pf10.3))')
914 & restyp(i),restyp(j),eps(i,j),sigma(i,j),chi(i,j),chi(j,i),
915 & chipp(i,j),chipp(j,i)
918 write (iout,'(/9a)') 'Residues',' alphasur1',' alphasur2',
919 & ' alphasur3',' alphasur4',' sigmap1',' sigmap2',
923 write (iout,'(2(a3,1x),8(0pf10.3))')
924 & restyp(i),restyp(j),(alphasur(k,i,j),k=1,4),
925 & sigmap1(i,j),sigmap2(j,i),chis(i,j),chis(j,i)
928 write (iout,'(/14a)') 'Residues',' nst ',' wst1',
929 & ' wst2',' wst3',' wst4',' dh11',' dh21',
930 & ' dh12',' dh22',' dt1',' dt2',' epsh1',
934 write (iout,'(2(a3,1x),i3,4f8.3,6f7.2,f9.5,f7.2)')
935 & restyp(i),restyp(j),nstate(i,j),(wstate(k,i,j),k=1,4),
936 & ((dhead(l,k,i,j),l=1,2),k=1,2),dtail(1,i,j),dtail(2,i,j),
937 & epshead(i,j),sig0head(i,j)
940 write (iout,'(/12a)') 'Residues',' ch1',' ch2',
941 & ' rborn1',' rborn2',' wqdip1',' wqdip2',
945 write (iout,'(2(a3,1x),2i4,5f10.3)')
946 & restyp(i),restyp(j),icharge(i),icharge(j),
947 & rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j)
950 write (iout,'(/12a)') 'Residues',
952 & ' alphpol2',' alphiso1',' alpiso2',
953 & ' alpiso3',' alpiso4',' sigiso1',' sigiso2',
957 write (iout,'(2(a3,1x),11f10.3)')
958 & restyp(i),restyp(j),alphapol(i,j),alphapol(j,i),
959 & (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(j,i),
967 C-----------------------------------------------------------------------
968 C Calculate the "working" parameters of SC interactions.
973 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
974 sigma(j,i)=sigma(i,j)
975 rs0(i,j)=dwa16*sigma(i,j)
982 write (iout,*) "IPOT=",ipot
983 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
984 & 'Working parameters of the SC interactions:',
985 & ' a ',' b ',' augm ',' sigma ',' r0 ',
990 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4 .or. ipot.eq.6 ) THEN
999 sigeps=dsign(1.0D0,epsij)
1001 aa(i,j)=epsij*rrij*rrij
1002 bb(i,j)=-sigeps*epsij*rrij
1005 IF ((ipot.gt.2).AND.(ipot.LT.6)) THEN
1006 sigt1sq=sigma0(i)**2
1007 sigt2sq=sigma0(j)**2
1010 ratsig1=sigt2sq/sigt1sq
1011 ratsig2=1.0D0/ratsig1
1012 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1013 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1014 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1018 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1019 sigmaii(i,j)=rsum_max
1020 sigmaii(j,i)=rsum_max
1022 c sigmaii(i,j)=r0(i,j)
1023 c sigmaii(j,i)=r0(i,j)
1025 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1026 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1027 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1028 augm(i,j)=epsij*r_augm**(2*expon)
1029 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1037 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1038 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1039 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1041 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3),2i3,10f8.4,
1043 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1044 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i),
1045 & icharge(i),icharge(j),chipp(i,j),chipp(j,i),
1046 & (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(j,i),
1047 & chis(i,j),chis(j,i),
1048 & nstate(i,j),(wstate(k,i,j),k=1,4),
1049 & ((dhead(l,k,i,j),l=1,2),k=1,2),dtail(1,i,j),dtail(2,i,j),
1050 & epshead(i,j),sig0head(i,j),
1051 & rborn(i,j),(wqdip(k,i,j),k=1,2),wquad(i,j),
1052 & alphapol(i,j),alphapol(j,i),
1053 & (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j)
1061 C Define the SC-p interaction constants
1065 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1067 c aad(i,1)=0.3D0*4.0D0**12
1068 C Following line for constants currently implemented
1069 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1070 aad(i,1)=1.5D0*4.0D0**12
1071 c aad(i,1)=0.17D0*5.6D0**12
1073 C "Soft" SC-p repulsion
1075 C Following line for constants currently implemented
1076 c aad(i,1)=0.3D0*4.0D0**6
1077 C "Hard" SC-p repulsion
1078 bad(i,1)=3.0D0*4.0D0**6
1079 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1088 C 8/9/01 Read the SC-p interaction constants from file
1091 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1094 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1095 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1096 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1097 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1101 write (iout,*) "Parameters of SC-p interactions:"
1103 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1104 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1109 C Define the constants of the disulfide bridge
1113 c Old arbitrary potential - commented out.
1118 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1119 c energy surface of diethyl disulfide.
1120 c A. Liwo and U. Kozlowska, 11/24/03
1131 write (iout,'(/a)') "Disulfide bridge parameters:"
1132 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1133 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1134 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1135 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1139 110 write (iout,*) "Error reading bond energy parameters."
1141 111 write (iout,*) "Error reading bending energy parameters."
1143 112 write (iout,*) "Error reading rotamer energy parameters."
1145 113 write (iout,*) "Error reading torsional energy parameters."
1147 114 write (iout,*) "Error reading double torsional energy parameters."
1150 & "Error reading cumulant (multibody energy) parameters."
1152 116 write (iout,*) "Error reading electrostatic energy parameters."
1154 117 write (iout,*) "Error reading side chain interaction parameters."
1156 118 write (iout,*) "Error reading SCp interaction parameters."
1158 119 write (iout,*) "Error reading SCCOR parameters"
1161 call MPI_Finalize(Ierror)