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,"SCCORPAR",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)
303 ithetyp(i)=-ithetyp(-i)
306 do i=-maxthetyp,maxthetyp
307 do j=-maxthetyp,maxthetyp
308 do k=-maxthetyp,maxthetyp
309 aa0thet(i,j,k,iblock)=0.0d0
311 aathet(l,i,j,k,iblock)=0.0d0
315 bbthet(m,l,i,j,k,iblock)=0.0d0
316 ccthet(m,l,i,j,k,iblock)=0.0d0
317 ddthet(m,l,i,j,k,iblock)=0.0d0
318 eethet(m,l,i,j,k,iblock)=0.0d0
324 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
325 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
335 do j=-nthetyp,nthetyp
336 do k=-nthetyp,nthetyp
337 read (ithep,'(6a)') res1
338 read (ithep,*) aa0thet(i,j,k,iblock)
339 read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
341 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
342 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
343 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
344 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle)
347 & (((ffthet(llll,lll,ll,i,j,k,iblock),
348 & ffthet(lll,llll,ll,i,j,k,iblock),
349 & ggthet(llll,lll,ll,i,j,k,iblock)
350 & ,ggthet(lll,llll,ll,i,j,k,iblock),
351 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
356 C For dummy ends assign glycine-type coefficients of theta-only terms; the
357 C coefficients of theta-and-gamma-dependent terms are zero.
362 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
363 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
365 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
366 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
369 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
371 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
374 C Substitution for D aminoacids from symmetry.
377 do j=-nthetyp,nthetyp
378 do k=-nthetyp,nthetyp
379 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
381 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
385 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
386 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
387 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
388 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
394 ffthet(llll,lll,ll,i,j,k,iblock)=
395 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
396 ffthet(lll,llll,ll,i,j,k,iblock)=
397 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
398 ggthet(llll,lll,ll,i,j,k,iblock)=
399 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
400 ggthet(lll,llll,ll,i,j,k,iblock)=
401 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
411 C Control printout of the coefficients of virtual-bond-angle potentials
414 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
418 write (iout,'(//4a)')
419 & 'Type ',onelett(i),onelett(j),onelett(k)
420 write (iout,'(//a,10x,a)') " l","a[l]"
421 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
422 write (iout,'(i2,1pe15.5)')
423 & (l,aathet(l,i,j,k),l=1,ntheterm)
425 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
426 & "b",l,"c",l,"d",l,"e",l
428 write (iout,'(i2,4(1pe15.5))') m,
429 & bbthet(m,l,i,j,k),ccthet(m,l,i,j,k,iblock),
430 & ddthet(m,l,i,j,k),eethet(m,l,i,j,k,iblock)
434 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
435 & "f+",l,"f-",l,"g+",l,"g-",l
438 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
439 & ffthet(n,m,l,i,j,k,iblock),
440 & ffthet(m,n,l,i,j,k,iblock),
441 & ggthet(n,m,l,i,j,k,iblock),
442 & ggthet(m,n,l,i,j,k,iblock)
455 C Read the parameters of the probability distribution/energy expression
456 C of the side chains.
459 cc write (iout,*) "tu dochodze",i
460 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
464 dsc_inv(i)=1.0D0/dsc(i)
475 read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
476 censc(1,1,-i)=censc(1,1,i)
477 censc(2,1,-i)=censc(2,1,i)
478 censc(3,1,-i)=-censc(3,1,i)
480 read (irotam,*) bsc(j,i)
481 read (irotam,*) (censc(k,j,i),k=1,3),
482 & ((blower(k,l,j),l=1,k),k=1,3)
483 censc(1,j,-i)=censc(1,j,i)
484 censc(2,j,-i)=censc(2,j,i)
485 censc(3,j,-i)=-censc(3,j,i)
486 C BSC is amplitude of Gaussian
493 akl=akl+blower(k,m,j)*blower(l,m,j)
497 if (((k.eq.3).and.(l.ne.3))
498 & .or.((l.eq.3).and.(k.ne.3))) then
499 gaussc(k,l,j,-i)=-akl
500 gaussc(l,k,j,-i)=-akl
512 write (iout,'(/a)') 'Parameters of side-chain local geometry'
516 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
517 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
518 c write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
519 c write (iout,'(a,f10.4,4(16x,f10.4))')
520 c & 'Center ',(bsc(j,i),j=1,nlobi)
521 c write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
522 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
523 & 'log h',(bsc(j,i),j=1,nlobi)
524 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
525 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
532 c blower(k,l,j)=gaussc(ind,j,i)
537 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
538 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
545 C Read scrot parameters for potentials determined from all-atom AM1 calculations
546 C added by Urszula Kozlowska 07/11/2007
554 read(irotam,*) sc_parmin(j,i)
562 C Read torsional parameters in old format
564 read (itorp,*) ntortyp,nterm_old
565 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
566 read (itorp,*) (itortyp(i),i=1,ntyp)
571 read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
577 write (iout,'(/a/)') 'Torsional constants:'
580 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
581 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
589 C Read torsional parameters
591 read (itorp,*) ntortyp
592 read (itorp,*) (itortyp(i),i=1,ntyp)
593 write (iout,*) 'ntortyp',ntortyp
596 itortyp(i)=-itortyp(-i)
598 c write (iout,*) 'ntortyp',ntortyp
600 do j=-ntortyp+1,ntortyp-1
601 read (itorp,*) nterm(i,j,iblock),
603 nterm(-i,-j,iblock)=nterm(i,j,iblock)
604 nlor(-i,-j,iblock)=nlor(i,j,iblock)
607 do k=1,nterm(i,j,iblock)
608 read (itorp,*) kk,v1(k,i,j,iblock),
610 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
611 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
612 v0ij=v0ij+si*v1(k,i,j,iblock)
615 do k=1,nlor(i,j,iblock)
616 read (itorp,*) kk,vlor1(k,i,j),
617 & vlor2(k,i,j),vlor3(k,i,j)
618 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
621 v0(-i,-j,iblock)=v0ij
627 write (iout,'(/a/)') 'Torsional constants:'
630 write (iout,*) 'ityp',i,' jtyp',j
631 write (iout,*) 'Fourier constants'
632 do k=1,nterm(i,j,iblock)
633 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
636 write (iout,*) 'Lorenz constants'
637 do k=1,nlor(i,j,iblock)
638 write (iout,'(3(1pe15.5))')
639 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
645 C 6/23/01 Read parameters for double torsionals
649 do j=-ntortyp+1,ntortyp-1
650 do k=-ntortyp+1,ntortyp-1
651 read (itordp,'(3a1)') t1,t2,t3
652 c write (iout,*) "OK onelett",
655 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
656 & .or. t3.ne.toronelet(k)) then
657 write (iout,*) "Error in double torsional parameter file",
660 call MPI_Finalize(Ierror)
662 stop "Error in double torsional parameter file"
664 read (itordp,*) ntermd_1(i,j,k,iblock),
665 & ntermd_2(i,j,k,iblock)
666 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
667 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
668 read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,
669 & ntermd_1(i,j,k,iblock))
670 read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,
671 & ntermd_1(i,j,k,iblock))
672 read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,
673 & ntermd_1(i,j,k,iblock))
674 read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,
675 & ntermd_1(i,j,k,iblock))
676 C Martix of D parameters for one dimesional foureir series
677 do l=1,ntermd_1(i,j,k,iblock)
678 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
679 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
680 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
681 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
682 c write(iout,*) "whcodze" ,
683 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
685 read (itordp,*) ((v2c(l,m,i,j,k,iblock),
686 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
687 & v2s(m,l,i,j,k,iblock),
688 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
689 C Martix of D parameters for two dimesional fourier series
690 do l=1,ntermd_2(i,j,k,iblock)
692 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
693 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
694 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
695 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
704 write (iout,*) 'Constants for double torsionals'
707 do j=-ntortyp+1,ntortyp-1
708 do k=-ntortyp+1,ntortyp-1
709 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
710 & ' nsingle',ntermd_1(i,j,k,iblock),
711 & ' ndouble',ntermd_2(i,j,k,iblock)
713 write (iout,*) 'Single angles:'
714 do l=1,ntermd_1(i,j,k,iblock)
715 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
716 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
717 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
718 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
721 write (iout,*) 'Pairs of angles:'
722 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
723 do l=1,ntermd_2(i,j,k,iblock)
724 write (iout,'(i5,20f10.5)')
725 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
728 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
729 do l=1,ntermd_2(i,j,k,iblock)
730 write (iout,'(i5,20f10.5)')
731 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
732 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
741 C Read of Side-chain backbone correlation parameters
742 C Modified 11 May 2012 by Adasko
745 read (isccor,*) nsccortyp
746 read (isccor,*) (isccortyp(i),i=1,ntyp)
748 isccortyp(i)=-isccortyp(-i)
750 iscprol=isccortyp(20)
751 c write (iout,*) 'ntortyp',ntortyp
753 cc maxinter is maximum interaction sites
758 &nterm_sccor(i,j),nlor_sccor(i,j)
759 write (iout,*) nterm_sccor(i,j)
765 nterm_sccor(-i,j)=nterm_sccor(i,j)
766 nterm_sccor(-i,-j)=nterm_sccor(i,j)
767 nterm_sccor(i,-j)=nterm_sccor(i,j)
768 write (iout,*) nterm_sccor(i,j),nterm_sccor(-i,j),
769 & nterm_sccor(-i,-j),nterm_sccor(i,-j)
770 do k=1,nterm_sccor(i,j)
771 read (isccor,*) kk,v1sccor(k,l,i,j)
773 if (j.eq.iscprol) then
774 if (i.eq.isccortyp(10)) then
775 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
776 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
778 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
779 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
780 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
781 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
782 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
783 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
784 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
785 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
788 if (i.eq.isccortyp(10)) then
789 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
790 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
792 if (j.eq.isccortyp(10)) then
793 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
794 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
796 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
797 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
798 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
799 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
800 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
801 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
805 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
806 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
807 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
808 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
811 do k=1,nlor_sccor(i,j)
812 read (isccor,*) kk,vlor1sccor(k,i,j),
813 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
814 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
815 &(1+vlor3sccor(k,i,j)**2)
817 v0sccor(l,i,j)=v0ijsccor
818 v0sccor(l,-i,j)=v0ijsccor1
819 v0sccor(l,i,-j)=v0ijsccor2
820 v0sccor(l,-i,-j)=v0ijsccor3
826 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
829 write (iout,*) 'ityp',i,' jtyp',j
830 write (iout,*) 'Fourier constants'
831 do k=1,nterm_sccor(i,j)
832 write (iout,'(2(1pe15.5))')
833 & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
835 write (iout,*) 'Lorenz constants'
836 do k=1,nlor_sccor(i,j)
837 write (iout,'(3(1pe15.5))')
838 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
844 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
845 C interaction energy of the Gly, Ala, and Pro prototypes.
847 read (ifourier,*) nloctyp
850 read (ifourier,*) (b(ii,i),ii=1,13)
852 write (iout,*) 'Type',i
853 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
861 B1tilde(1,i) = b(3,i)
862 B1tilde(2,i) =-b(5,i)
863 B1tilde(1,-i) =-b(3,i)
864 B1tilde(2,-i) =b(5,i)
888 Ctilde(2,1,i)=-b(9,i)
890 Ctilde(1,1,-i)=b(7,i)
891 Ctilde(1,2,-i)=-b(9,i)
892 Ctilde(2,1,-i)=b(9,i)
893 Ctilde(2,2,-i)=b(7,i)
895 c Ctilde(1,1,i)=0.0d0
896 c Ctilde(1,2,i)=0.0d0
897 c Ctilde(2,1,i)=0.0d0
898 c Ctilde(2,2,i)=0.0d0
913 Dtilde(2,1,i)=-b(8,i)
915 Dtilde(1,1,-i)=b(6,i)
916 Dtilde(1,2,-i)=-b(8,i)
917 Dtilde(2,1,-i)=b(8,i)
918 Dtilde(2,2,-i)=b(6,i)
920 c Dtilde(1,1,i)=0.0d0
921 c Dtilde(1,2,i)=0.0d0
922 c Dtilde(2,1,i)=0.0d0
923 c Dtilde(2,2,i)=0.0d0
924 EE(1,1,i)= b(10,i)+b(11,i)
925 EE(2,2,i)=-b(10,i)+b(11,i)
926 EE(2,1,i)= b(12,i)-b(13,i)
927 EE(1,2,i)= b(12,i)+b(13,i)
928 EE(1,1,-i)= b(10,i)+b(11,i)
929 EE(2,2,-i)=-b(10,i)+b(11,i)
930 EE(2,1,-i)=-b(12,i)+b(13,i)
931 EE(1,2,-i)=-b(12,i)-b(13,i)
937 c ee(2,1,i)=ee(1,2,i)
942 write (iout,*) 'Type',i
944 c write (iout,'(f10.5)') B1(:,i)
945 write(iout,*) B1(1,i),B1(2,i)
947 c write (iout,'(f10.5)') B2(:,i)
948 write(iout,*) B2(1,i),B2(2,i)
951 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
955 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
959 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
964 C Read electrostatic-interaction parameters
967 write (iout,'(/a)') 'Electrostatic interaction constants:'
968 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
969 & 'IT','JT','APP','BPP','AEL6','AEL3'
971 read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
972 read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
973 read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
974 read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
979 app (i,j)=epp(i,j)*rri*rri
980 bpp (i,j)=-2.0D0*epp(i,j)*rri
981 ael6(i,j)=elpp6(i,j)*4.2D0**6
982 ael3(i,j)=elpp3(i,j)*4.2D0**3
983 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
984 & ael6(i,j),ael3(i,j)
988 C Read side-chain interaction parameters.
990 read (isidep,*) ipot,expon
991 if (ipot.lt.1 .or. ipot.gt.5) then
992 write (iout,'(2a)') 'Error while reading SC interaction',
993 & 'potential file - unknown potential type.'
997 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
998 & ', exponents are ',expon,2*expon
999 goto (10,20,30,30,40) ipot
1000 C----------------------- LJ potential ---------------------------------
1001 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1003 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1004 write (iout,'(a/)') 'The epsilon array:'
1005 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1006 write (iout,'(/a)') 'One-body parameters:'
1007 write (iout,'(a,4x,a)') 'residue','sigma'
1008 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1011 C----------------------- LJK potential --------------------------------
1012 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1013 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1015 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1016 write (iout,'(a/)') 'The epsilon array:'
1017 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1018 write (iout,'(/a)') 'One-body parameters:'
1019 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1020 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1024 C---------------------- GB or BP potential -----------------------------
1025 30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1026 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
1028 C For the GB potential convert sigma'**2 into chi'
1031 chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
1035 write (iout,'(/a/)') 'Parameters of the BP potential:'
1036 write (iout,'(a/)') 'The epsilon array:'
1037 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1038 write (iout,'(/a)') 'One-body parameters:'
1039 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1041 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1042 & chip(i),alp(i),i=1,ntyp)
1045 C--------------------- GBV potential -----------------------------------
1046 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1047 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1048 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1050 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1051 write (iout,'(a/)') 'The epsilon array:'
1052 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1053 write (iout,'(/a)') 'One-body parameters:'
1054 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1055 & 's||/s_|_^2',' chip ',' alph '
1056 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1057 & sigii(i),chip(i),alp(i),i=1,ntyp)
1061 C-----------------------------------------------------------------------
1062 C Calculate the "working" parameters of SC interactions.
1070 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1071 sigma(j,i)=sigma(i,j)
1072 rs0(i,j)=dwa16*sigma(i,j)
1076 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1077 & 'Working parameters of the SC interactions:',
1078 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1083 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1092 sigeps=dsign(1.0D0,epsij)
1094 aa(i,j)=epsij*rrij*rrij
1095 bb(i,j)=-sigeps*epsij*rrij
1099 sigt1sq=sigma0(i)**2
1100 sigt2sq=sigma0(j)**2
1103 ratsig1=sigt2sq/sigt1sq
1104 ratsig2=1.0D0/ratsig1
1105 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1106 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1107 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1111 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1112 sigmaii(i,j)=rsum_max
1113 sigmaii(j,i)=rsum_max
1115 c sigmaii(i,j)=r0(i,j)
1116 c sigmaii(j,i)=r0(i,j)
1118 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1119 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1120 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1121 augm(i,j)=epsij*r_augm**(2*expon)
1122 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1129 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1130 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1131 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1136 C Define the SC-p interaction constants
1140 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1142 c aad(i,1)=0.3D0*4.0D0**12
1143 C Following line for constants currently implemented
1144 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1145 aad(i,1)=1.5D0*4.0D0**12
1146 c aad(i,1)=0.17D0*5.6D0**12
1148 C "Soft" SC-p repulsion
1150 C Following line for constants currently implemented
1151 c aad(i,1)=0.3D0*4.0D0**6
1152 C "Hard" SC-p repulsion
1153 bad(i,1)=3.0D0*4.0D0**6
1154 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1163 C 8/9/01 Read the SC-p interaction constants from file
1166 read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1169 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1170 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1171 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1172 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1176 write (iout,*) "Parameters of SC-p interactions:"
1178 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1179 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1184 C Define the constants of the disulfide bridge
1188 c Old arbitrary potential - commented out.
1193 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1194 c energy surface of diethyl disulfide.
1195 c A. Liwo and U. Kozlowska, 11/24/03
1206 write (iout,'(/a)') "Disulfide bridge parameters:"
1207 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1208 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1209 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1210 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,