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"/
26 character*1 toronelet(-2:2)/"p","a","G","A","P"/
28 dimension blower(3,3,maxlob)
29 character*800 controlcard
30 character*256 bondname_t,thetname_t,rotname_t,torname_t,
31 & tordname_t,fouriername_t,elename_t,sidename_t,scpname_t,
37 double precision ip,mp
41 C Set LPRINT=.TRUE. for debugging
42 dwa16=2.0d0**(1.0d0/6.0d0)
45 C Assign virtual-bond length
49 call card_concat(controlcard,.true.)
52 key = wname(i)(:ilen(wname(i)))
53 call reada(controlcard,key(:ilen(key)),ww(i),1.0d0)
56 write (iout,*) "iparm",iparm," myparm",myparm
57 c If reading not own parameters, skip assignment
59 if (iparm.eq.myparm .or. .not.separate_parset) then
62 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,*) vbldp0,akp
156 read (ibond,*) vbldsc0(1,i),aksc(1,i)
157 dsc(i) = vbldsc0(1,i)
161 dsc_inv(i)=1.0D0/dsc(i)
165 read (ibond,*) ijunk,vbldp0,akp,rjunk
167 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
169 dsc(i) = vbldsc0(1,i)
173 dsc_inv(i)=1.0D0/dsc(i)
178 write(iout,'(/a/)')"Force constants virtual bonds:"
179 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
181 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
183 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
184 & vbldsc0(1,i),aksc(1,i),abond0(1,i)
186 write (iout,'(13x,3f10.5)')
187 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
193 C Read the parameters of the probability distribution/energy expression
194 C of the virtual-bond valence angles theta
197 read (ithep,*) a0thet(i),(athet(j,i,1,1),j=1,2),
198 & (bthet(j,i,1,1),j=1,2)
199 read (ithep,*) (polthet(j,i),j=0,3)
200 read (ithep,*) (gthet(j,i),j=1,3)
201 read (ithep,*) theta0(i),sig0(i),sigc0(i)
205 athet(1,i,1,-1)=athet(1,i,1,1)
206 athet(2,i,1,-1)=athet(2,i,1,1)
207 bthet(1,i,1,-1)=-bthet(1,i,1,1)
208 bthet(2,i,1,-1)=-bthet(2,i,1,1)
209 athet(1,i,-1,1)=-athet(1,i,1,1)
210 athet(2,i,-1,1)=-athet(2,i,1,1)
211 bthet(1,i,-1,1)=bthet(1,i,1,1)
212 bthet(2,i,-1,1)=bthet(2,i,1,1)
216 athet(1,i,-1,-1)=athet(1,-i,1,1)
217 athet(2,i,-1,-1)=-athet(2,-i,1,1)
218 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
219 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
220 athet(1,i,-1,1)=athet(1,-i,1,1)
221 athet(2,i,-1,1)=-athet(2,-i,1,1)
222 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
223 bthet(2,i,-1,1)=bthet(2,-i,1,1)
224 athet(1,i,1,-1)=-athet(1,-i,1,1)
225 athet(2,i,1,-1)=athet(2,-i,1,1)
226 bthet(1,i,1,-1)=bthet(1,-i,1,1)
227 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
232 polthet(j,i)=polthet(j,-i)
235 gthet(j,i)=gthet(j,-i)
241 c & 'Parameters of the virtual-bond valence angles:'
242 c write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
243 c & ' ATHETA0 ',' A1 ',' A2 ',
246 c write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
247 c & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
249 c write (iout,'(/a/9x,5a/79(1h-))')
250 c & 'Parameters of the expression for sigma(theta_c):',
251 c & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
252 c & ' ALPH3 ',' SIGMA0C '
254 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
255 c & (polthet(j,i),j=0,3),sigc0(i)
257 c write (iout,'(/a/9x,5a/79(1h-))')
258 c & 'Parameters of the second gaussian:',
259 c & ' THETA0 ',' SIGMA0 ',' G1 ',
262 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
263 c & sig0(i),(gthet(j,i),j=1,3)
266 & 'Parameters of the virtual-bond valence angles:'
267 write (iout,'(/a/9x,5a/79(1h-))')
268 & 'Coefficients of expansion',
269 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
270 & ' b1*10^1 ',' b2*10^1 '
272 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
273 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
274 & (10*bthet(j,i,1,1),j=1,2)
276 write (iout,'(/a/9x,5a/79(1h-))')
277 & 'Parameters of the expression for sigma(theta_c):',
278 & ' alpha0 ',' alph1 ',' alph2 ',
279 & ' alhp3 ',' sigma0c '
281 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
282 & (polthet(j,i),j=0,3),sigc0(i)
284 write (iout,'(/a/9x,5a/79(1h-))')
285 & 'Parameters of the second gaussian:',
286 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
289 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
290 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
295 C Read the parameters of Utheta determined from ab initio surfaces
296 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
298 read (ithep,*) nthetyp,ntheterm,ntheterm2,
299 & ntheterm3,nsingle,ndouble
300 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
301 read (ithep,*) (ithetyp(i),i=1,ntyp1)
307 aathet(l,i,j,k)=0.0d0
311 bbthet(m,l,i,j,k)=0.0d0
312 ccthet(m,l,i,j,k)=0.0d0
313 ddthet(m,l,i,j,k)=0.0d0
314 eethet(m,l,i,j,k)=0.0d0
320 ffthet(mm,m,l,i,j,k)=0.0d0
321 ggthet(mm,m,l,i,j,k)=0.0d0
331 read (ithep,'(3a)') res1,res2,res3
332 read (ithep,*) aa0thet(i,j,k)
333 read (ithep,*)(aathet(l,i,j,k),l=1,ntheterm)
335 & ((bbthet(lll,ll,i,j,k),lll=1,nsingle),
336 & (ccthet(lll,ll,i,j,k),lll=1,nsingle),
337 & (ddthet(lll,ll,i,j,k),lll=1,nsingle),
338 & (eethet(lll,ll,i,j,k),lll=1,nsingle),ll=1,ntheterm2)
340 & (((ffthet(llll,lll,ll,i,j,k),ffthet(lll,llll,ll,i,j,k),
341 & ggthet(llll,lll,ll,i,j,k),ggthet(lll,llll,ll,i,j,k),
342 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
347 C For dummy ends assign glycine-type coefficients of theta-only terms; the
348 C coefficients of theta-and-gamma-dependent terms are zero.
353 aathet(l,i,j,nthetyp+1)=aathet(l,i,j,1)
354 aathet(l,nthetyp+1,i,j)=aathet(l,1,i,j)
356 aa0thet(i,j,nthetyp+1)=aa0thet(i,j,1)
357 aa0thet(nthetyp+1,i,j)=aa0thet(1,i,j)
360 aathet(l,nthetyp+1,i,nthetyp+1)=aathet(l,1,i,1)
362 aa0thet(nthetyp+1,i,nthetyp+1)=aa0thet(1,i,1)
365 C Control printout of the coefficients of virtual-bond-angle potentials
368 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
372 write (iout,'(//4a)')
373 & 'Type ',onelett(i),onelett(j),onelett(k)
374 write (iout,'(//a,10x,a)') " l","a[l]"
375 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k)
376 write (iout,'(i2,1pe15.5)')
377 & (l,aathet(l,i,j,k),l=1,ntheterm)
379 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
380 & "b",l,"c",l,"d",l,"e",l
382 write (iout,'(i2,4(1pe15.5))') m,
383 & bbthet(m,l,i,j,k),ccthet(m,l,i,j,k),
384 & ddthet(m,l,i,j,k),eethet(m,l,i,j,k)
388 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
389 & "f+",l,"f-",l,"g+",l,"g-",l
392 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
393 & ffthet(n,m,l,i,j,k),ffthet(m,n,l,i,j,k),
394 & ggthet(n,m,l,i,j,k),ggthet(m,n,l,i,j,k)
407 C Read the parameters of the probability distribution/energy expression
408 C of the side chains.
411 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
415 dsc_inv(i)=1.0D0/dsc(i)
426 read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
427 censc(1,1,-i)=censc(1,1,i)
428 censc(2,1,-i)=censc(2,1,i)
429 censc(3,1,-i)=-censc(3,1,i)
431 read (irotam,*) bsc(j,i)
432 read (irotam,*) (censc(k,j,i),k=1,3),
433 & ((blower(k,l,j),l=1,k),k=1,3)
434 censc(1,j,-i)=censc(1,j,i)
435 censc(2,j,-i)=censc(2,j,i)
436 censc(3,j,-i)=-censc(3,j,i)
443 akl=akl+blower(k,m,j)*blower(l,m,j)
447 if (((k.eq.3).and.(l.ne.3))
448 & .or.((l.eq.3).and.(k.ne.3))) then
449 gaussc(k,l,j,-i)=-akl
450 gaussc(l,k,j,-i)=-akl
462 write (iout,'(/a)') 'Parameters of side-chain local geometry'
466 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
467 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
468 c write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
469 c write (iout,'(a,f10.4,4(16x,f10.4))')
470 c & 'Center ',(bsc(j,i),j=1,nlobi)
471 c write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
472 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
473 & 'log h',(bsc(j,i),j=1,nlobi)
474 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
475 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
482 c blower(k,l,j)=gaussc(ind,j,i)
487 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
488 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
495 C Read scrot parameters for potentials determined from all-atom AM1 calculations
496 C added by Urszula Kozlowska 07/11/2007
504 read(irotam,*) sc_parmin(j,i)
512 C Read torsional parameters in old format
514 read (itorp,*) ntortyp,nterm_old
515 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
516 read (itorp,*) (itortyp(i),i=1,ntyp)
521 read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
527 write (iout,'(/a/)') 'Torsional constants:'
530 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
531 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
539 C Read torsional parameters
541 read (itorp,*) ntortyp
542 read (itorp,*) (itortyp(i),i=1,ntyp)
543 write (iout,*) 'ntortyp',ntortyp
546 read (itorp,*) nterm(i,j),nlor(i,j)
550 read (itorp,*) kk,v1(k,i,j),v2(k,i,j)
551 v0ij=v0ij+si*v1(k,i,j)
555 read (itorp,*) kk,vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
556 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
563 write (iout,'(/a/)') 'Torsional constants:'
566 write (iout,*) 'ityp',i,' jtyp',j
567 write (iout,*) 'Fourier constants'
569 write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
571 write (iout,*) 'Lorenz constants'
573 write (iout,'(3(1pe15.5))')
574 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
580 C 6/23/01 Read parameters for double torsionals
583 do j=-ntortyp,ntortyp
585 read (itordp,'(3a1)') t1,t2,t3
586 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
587 & .or. t3.ne.toronelet(k)) then
588 write (iout,*) "Error in double torsional parameter file",
590 stop "Error in double torsional parameter file"
592 read (itordp,*) ntermd_1(i,j,k),ntermd_2(i,j,k)
593 read (itordp,*) (v1c(1,l,i,j,k),l=1,ntermd_1(i,j,k))
594 read (itordp,*) (v1s(1,l,i,j,k),l=1,ntermd_1(i,j,k))
595 read (itordp,*) (v1c(2,l,i,j,k),l=1,ntermd_1(i,j,k))
596 read (itordp,*) (v1s(2,l,i,j,k),l=1,ntermd_1(i,j,k))
597 read (itordp,*) ((v2c(l,m,i,j,k),v2c(m,l,i,j,k),
598 & v2s(l,m,i,j,k),v2s(m,l,i,j,k),m=1,l-1),l=1,ntermd_2(i,j,k))
604 write (iout,*) 'Constants for double torsionals'
608 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
609 & ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k)
611 write (iout,*) 'Single angles:'
612 do l=1,ntermd_1(i,j,k)
613 write (iout,'(i5,2f10.5,5x,2f10.5)') l,
614 & v1c(1,l,i,j,k),v1s(1,l,i,j,k),
615 & v1c(2,l,i,j,k),v1s(2,l,i,j,k)
618 write (iout,*) 'Pairs of angles:'
619 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
620 do l=1,ntermd_2(i,j,k)
621 write (iout,'(i5,20f10.5)')
622 & l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k))
625 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
626 do l=1,ntermd_2(i,j,k)
627 write (iout,'(i5,20f10.5)')
628 & l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k))
636 C Read of Side-chain backbone correlation parameters
637 C Modified 11 May 2012 by Adasko
640 read (isccor,*) nsccortyp
641 read (isccor,*) (isccortyp(i),i=1,ntyp)
642 c write (iout,*) 'ntortyp',ntortyp
644 cc maxinter is maximum interaction sites
648 read (isccor,*) nterm_sccor(i,j),nlor_sccor(i,j)
652 do k=1,nterm_sccor(i,j)
653 read (isccor,*) kk,v1sccor(k,l,i,j)
655 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
658 do k=1,nlor_sccor(i,j)
659 read (isccor,*) kk,vlor1sccor(k,i,j),
660 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
661 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
662 &(1+vlor3sccor(k,i,j)**2)
664 v0sccor(i,j)=v0ijsccor
671 write (iout,'(/a/)') 'Torsional constants:'
674 write (iout,*) 'ityp',i,' jtyp',j
675 write (iout,*) 'Fourier constants'
676 do k=1,nterm_sccor(i,j)
677 write (iout,'(2(1pe15.5))')v1sccor(k,l,i,j),v2sccor(k,l,i,j)
679 write (iout,*) 'Lorenz constants'
680 do k=1,nlor_sccor(i,j)
681 write (iout,'(3(1pe15.5))')
682 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
689 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
690 C interaction energy of the Gly, Ala, and Pro prototypes.
692 read (ifourier,*) nloctyp
694 read (ifourier,*,end=115,err=115)
695 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
697 write (iout,*) 'Type',i
698 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
740 c Ctilde(1,1,i)=0.0d0
741 c Ctilde(1,2,i)=0.0d0
742 c Ctilde(2,1,i)=0.0d0
743 c Ctilde(2,2,i)=0.0d0
765 c Dtilde(1,1,i)=0.0d0
766 c Dtilde(1,2,i)=0.0d0
767 c Dtilde(2,1,i)=0.0d0
768 c Dtilde(2,2,i)=0.0d0
769 EE(1,1,i)= b(10)+b(11)
770 EE(2,2,i)=-b(10)+b(11)
771 EE(2,1,i)= b(12)-b(13)
772 EE(1,2,i)= b(12)+b(13)
773 EE(1,1,-i)= b(10)+b(11)
774 EE(2,2,-i)=-b(10)+b(11)
775 EE(2,1,-i)=-b(12)+b(13)
776 EE(1,2,-i)=-b(12)-b(13)
782 c ee(2,1,i)=ee(1,2,i)
786 write (iout,*) 'Type',i
788 c write (iout,'(f10.5)') B1(:,i)
789 write(iout,*) B1(1,i),B1(2,i)
791 c write (iout,'(f10.5)') B2(:,i)
792 write(iout,*) B2(1,i),B2(2,i)
795 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
799 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
803 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
808 C Read electrostatic-interaction parameters
811 write (iout,'(/a)') 'Electrostatic interaction constants:'
812 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
813 & 'IT','JT','APP','BPP','AEL6','AEL3'
815 read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
816 read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
817 read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
818 read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
823 app (i,j)=epp(i,j)*rri*rri
824 bpp (i,j)=-2.0D0*epp(i,j)*rri
825 ael6(i,j)=elpp6(i,j)*4.2D0**6
826 ael3(i,j)=elpp3(i,j)*4.2D0**3
827 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
828 & ael6(i,j),ael3(i,j)
832 C Read side-chain interaction parameters.
834 read (isidep,*) ipot,expon
835 if (ipot.lt.1 .or. ipot.gt.5) then
836 write (iout,'(2a)') 'Error while reading SC interaction',
837 & 'potential file - unknown potential type.'
841 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
842 & ', exponents are ',expon,2*expon
843 goto (10,20,30,30,40) ipot
844 C----------------------- LJ potential ---------------------------------
845 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
847 write (iout,'(/a/)') 'Parameters of the LJ potential:'
848 write (iout,'(a/)') 'The epsilon array:'
849 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
850 write (iout,'(/a)') 'One-body parameters:'
851 write (iout,'(a,4x,a)') 'residue','sigma'
852 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
855 C----------------------- LJK potential --------------------------------
856 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
857 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
859 write (iout,'(/a/)') 'Parameters of the LJK potential:'
860 write (iout,'(a/)') 'The epsilon array:'
861 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
862 write (iout,'(/a)') 'One-body parameters:'
863 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
864 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
868 C---------------------- GB or BP potential -----------------------------
869 30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
870 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
872 C For the GB potential convert sigma'**2 into chi'
875 chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
879 write (iout,'(/a/)') 'Parameters of the BP potential:'
880 write (iout,'(a/)') 'The epsilon array:'
881 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
882 write (iout,'(/a)') 'One-body parameters:'
883 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
885 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
886 & chip(i),alp(i),i=1,ntyp)
889 C--------------------- GBV potential -----------------------------------
890 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
891 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
892 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
894 write (iout,'(/a/)') 'Parameters of the GBV potential:'
895 write (iout,'(a/)') 'The epsilon array:'
896 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
897 write (iout,'(/a)') 'One-body parameters:'
898 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
899 & 's||/s_|_^2',' chip ',' alph '
900 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
901 & sigii(i),chip(i),alp(i),i=1,ntyp)
905 C-----------------------------------------------------------------------
906 C Calculate the "working" parameters of SC interactions.
914 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
915 sigma(j,i)=sigma(i,j)
916 rs0(i,j)=dwa16*sigma(i,j)
920 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
921 & 'Working parameters of the SC interactions:',
922 & ' a ',' b ',' augm ',' sigma ',' r0 ',
927 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
936 sigeps=dsign(1.0D0,epsij)
938 aa(i,j)=epsij*rrij*rrij
939 bb(i,j)=-sigeps*epsij*rrij
947 ratsig1=sigt2sq/sigt1sq
948 ratsig2=1.0D0/ratsig1
949 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
950 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
951 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
955 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
956 sigmaii(i,j)=rsum_max
957 sigmaii(j,i)=rsum_max
959 c sigmaii(i,j)=r0(i,j)
960 c sigmaii(j,i)=r0(i,j)
962 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
963 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
964 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
965 augm(i,j)=epsij*r_augm**(2*expon)
966 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
973 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
974 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
975 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
980 C Define the SC-p interaction constants
984 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
986 c aad(i,1)=0.3D0*4.0D0**12
987 C Following line for constants currently implemented
988 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
989 aad(i,1)=1.5D0*4.0D0**12
990 c aad(i,1)=0.17D0*5.6D0**12
992 C "Soft" SC-p repulsion
994 C Following line for constants currently implemented
995 c aad(i,1)=0.3D0*4.0D0**6
996 C "Hard" SC-p repulsion
997 bad(i,1)=3.0D0*4.0D0**6
998 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1007 C 8/9/01 Read the SC-p interaction constants from file
1010 read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1013 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1014 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1015 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1016 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1020 write (iout,*) "Parameters of SC-p interactions:"
1022 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1023 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1028 C Define the constants of the disulfide bridge
1032 c Old arbitrary potential - commented out.
1037 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1038 c energy surface of diethyl disulfide.
1039 c A. Liwo and U. Kozlowska, 11/24/03
1050 write (iout,'(/a)') "Disulfide bridge parameters:"
1051 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1052 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1053 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1054 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,