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,vbldpdum,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,vbldpdum,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)
191 read(iliptranpar,*) pepliptran
193 read(iliptranpar,*) liptranene(i)
198 C Read the parameters of the probability distribution/energy expression
199 C of the virtual-bond valence angles theta
202 read (ithep,*) a0thet(i),(athet(j,i,1,1),j=1,2),
203 & (bthet(j,i,1,1),j=1,2)
204 read (ithep,*) (polthet(j,i),j=0,3)
205 read (ithep,*) (gthet(j,i),j=1,3)
206 read (ithep,*) theta0(i),sig0(i),sigc0(i)
210 athet(1,i,1,-1)=athet(1,i,1,1)
211 athet(2,i,1,-1)=athet(2,i,1,1)
212 bthet(1,i,1,-1)=-bthet(1,i,1,1)
213 bthet(2,i,1,-1)=-bthet(2,i,1,1)
214 athet(1,i,-1,1)=-athet(1,i,1,1)
215 athet(2,i,-1,1)=-athet(2,i,1,1)
216 bthet(1,i,-1,1)=bthet(1,i,1,1)
217 bthet(2,i,-1,1)=bthet(2,i,1,1)
221 athet(1,i,-1,-1)=athet(1,-i,1,1)
222 athet(2,i,-1,-1)=-athet(2,-i,1,1)
223 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
224 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
225 athet(1,i,-1,1)=athet(1,-i,1,1)
226 athet(2,i,-1,1)=-athet(2,-i,1,1)
227 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
228 bthet(2,i,-1,1)=bthet(2,-i,1,1)
229 athet(1,i,1,-1)=-athet(1,-i,1,1)
230 athet(2,i,1,-1)=athet(2,-i,1,1)
231 bthet(1,i,1,-1)=bthet(1,-i,1,1)
232 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
237 polthet(j,i)=polthet(j,-i)
240 gthet(j,i)=gthet(j,-i)
246 c & 'Parameters of the virtual-bond valence angles:'
247 c write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
248 c & ' ATHETA0 ',' A1 ',' A2 ',
251 c write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
252 c & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
254 c write (iout,'(/a/9x,5a/79(1h-))')
255 c & 'Parameters of the expression for sigma(theta_c):',
256 c & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
257 c & ' ALPH3 ',' SIGMA0C '
259 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
260 c & (polthet(j,i),j=0,3),sigc0(i)
262 c write (iout,'(/a/9x,5a/79(1h-))')
263 c & 'Parameters of the second gaussian:',
264 c & ' THETA0 ',' SIGMA0 ',' G1 ',
267 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
268 c & sig0(i),(gthet(j,i),j=1,3)
271 & 'Parameters of the virtual-bond valence angles:'
272 write (iout,'(/a/9x,5a/79(1h-))')
273 & 'Coefficients of expansion',
274 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
275 & ' b1*10^1 ',' b2*10^1 '
277 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
278 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
279 & (10*bthet(j,i,1,1),j=1,2)
281 write (iout,'(/a/9x,5a/79(1h-))')
282 & 'Parameters of the expression for sigma(theta_c):',
283 & ' alpha0 ',' alph1 ',' alph2 ',
284 & ' alhp3 ',' sigma0c '
286 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
287 & (polthet(j,i),j=0,3),sigc0(i)
289 write (iout,'(/a/9x,5a/79(1h-))')
290 & 'Parameters of the second gaussian:',
291 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
294 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
295 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
300 C Read the parameters of Utheta determined from ab initio surfaces
301 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
303 c write (iout,*) "tu dochodze"
304 read (ithep,*) nthetyp,ntheterm,ntheterm2,
305 & ntheterm3,nsingle,ndouble
306 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
307 read (ithep,*) (ithetyp(i),i=1,ntyp1)
309 ithetyp(i)=-ithetyp(-i)
311 c write (iout,*) "tu dochodze"
313 do i=-maxthetyp,maxthetyp
314 do j=-maxthetyp,maxthetyp
315 do k=-maxthetyp,maxthetyp
316 aa0thet(i,j,k,iblock)=0.0d0
318 aathet(l,i,j,k,iblock)=0.0d0
322 bbthet(m,l,i,j,k,iblock)=0.0d0
323 ccthet(m,l,i,j,k,iblock)=0.0d0
324 ddthet(m,l,i,j,k,iblock)=0.0d0
325 eethet(m,l,i,j,k,iblock)=0.0d0
331 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
332 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
342 do j=-nthetyp,nthetyp
343 do k=-nthetyp,nthetyp
344 read (ithep,'(6a)') res1
345 read (ithep,*) aa0thet(i,j,k,iblock)
346 read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
348 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
349 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
350 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
351 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle)
354 & (((ffthet(llll,lll,ll,i,j,k,iblock),
355 & ffthet(lll,llll,ll,i,j,k,iblock),
356 & ggthet(llll,lll,ll,i,j,k,iblock)
357 & ,ggthet(lll,llll,ll,i,j,k,iblock),
358 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
363 C For dummy ends assign glycine-type coefficients of theta-only terms; the
364 C coefficients of theta-and-gamma-dependent terms are zero.
369 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
370 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
372 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
373 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
376 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
378 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
381 C Substitution for D aminoacids from symmetry.
384 do j=-nthetyp,nthetyp
385 do k=-nthetyp,nthetyp
386 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
388 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
392 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
393 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
394 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
395 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
401 ffthet(llll,lll,ll,i,j,k,iblock)=
402 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
403 ffthet(lll,llll,ll,i,j,k,iblock)=
404 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
405 ggthet(llll,lll,ll,i,j,k,iblock)=
406 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
407 ggthet(lll,llll,ll,i,j,k,iblock)=
408 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
418 C Control printout of the coefficients of virtual-bond-angle potentials
421 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
425 write (iout,'(//4a)')
426 & 'Type ',onelett(i),onelett(j),onelett(k)
427 write (iout,'(//a,10x,a)') " l","a[l]"
428 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
429 write (iout,'(i2,1pe15.5)')
430 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
432 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
433 & "b",l,"c",l,"d",l,"e",l
435 write (iout,'(i2,4(1pe15.5))') m,
436 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
437 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
441 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
442 & "f+",l,"f-",l,"g+",l,"g-",l
445 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
446 & ffthet(n,m,l,i,j,k,iblock),
447 & ffthet(m,n,l,i,j,k,iblock),
448 & ggthet(n,m,l,i,j,k,iblock),
449 & ggthet(m,n,l,i,j,k,iblock)
462 C Read the parameters of the probability distribution/energy expression
463 C of the side chains.
466 cc write (iout,*) "tu dochodze",i
467 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
471 dsc_inv(i)=1.0D0/dsc(i)
482 read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
483 censc(1,1,-i)=censc(1,1,i)
484 censc(2,1,-i)=censc(2,1,i)
485 censc(3,1,-i)=-censc(3,1,i)
487 read (irotam,*) bsc(j,i)
488 read (irotam,*) (censc(k,j,i),k=1,3),
489 & ((blower(k,l,j),l=1,k),k=1,3)
490 censc(1,j,-i)=censc(1,j,i)
491 censc(2,j,-i)=censc(2,j,i)
492 censc(3,j,-i)=-censc(3,j,i)
493 C BSC is amplitude of Gaussian
500 akl=akl+blower(k,m,j)*blower(l,m,j)
504 if (((k.eq.3).and.(l.ne.3))
505 & .or.((l.eq.3).and.(k.ne.3))) then
506 gaussc(k,l,j,-i)=-akl
507 gaussc(l,k,j,-i)=-akl
519 write (iout,'(/a)') 'Parameters of side-chain local geometry'
523 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
524 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
525 c write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
526 c write (iout,'(a,f10.4,4(16x,f10.4))')
527 c & 'Center ',(bsc(j,i),j=1,nlobi)
528 c write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
529 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
530 & 'log h',(bsc(j,i),j=1,nlobi)
531 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
532 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
539 c blower(k,l,j)=gaussc(ind,j,i)
544 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
545 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
552 C Read scrot parameters for potentials determined from all-atom AM1 calculations
553 C added by Urszula Kozlowska 07/11/2007
561 read(irotam,*) sc_parmin(j,i)
569 C Read torsional parameters in old format
571 read (itorp,*) ntortyp,nterm_old
572 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
573 read (itorp,*) (itortyp(i),i=1,ntyp)
578 read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
584 write (iout,'(/a/)') 'Torsional constants:'
587 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
588 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
596 C Read torsional parameters
598 read (itorp,*) ntortyp
599 read (itorp,*) (itortyp(i),i=1,ntyp)
600 write (iout,*) 'ntortyp',ntortyp
603 itortyp(i)=-itortyp(-i)
605 c write (iout,*) 'ntortyp',ntortyp
607 do j=-ntortyp+1,ntortyp-1
608 read (itorp,*) nterm(i,j,iblock),
610 nterm(-i,-j,iblock)=nterm(i,j,iblock)
611 nlor(-i,-j,iblock)=nlor(i,j,iblock)
614 do k=1,nterm(i,j,iblock)
615 read (itorp,*) kk,v1(k,i,j,iblock),
617 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
618 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
619 v0ij=v0ij+si*v1(k,i,j,iblock)
622 do k=1,nlor(i,j,iblock)
623 read (itorp,*) kk,vlor1(k,i,j),
624 & vlor2(k,i,j),vlor3(k,i,j)
625 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
628 v0(-i,-j,iblock)=v0ij
634 write (iout,'(/a/)') 'Torsional constants:'
637 write (iout,*) 'ityp',i,' jtyp',j
638 write (iout,*) 'Fourier constants'
639 do k=1,nterm(i,j,iblock)
640 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
643 write (iout,*) 'Lorenz constants'
644 do k=1,nlor(i,j,iblock)
645 write (iout,'(3(1pe15.5))')
646 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
652 C 6/23/01 Read parameters for double torsionals
656 do j=-ntortyp+1,ntortyp-1
657 do k=-ntortyp+1,ntortyp-1
658 read (itordp,'(3a1)') t1,t2,t3
659 c write (iout,*) "OK onelett",
662 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
663 & .or. t3.ne.toronelet(k)) then
664 write (iout,*) "Error in double torsional parameter file",
667 call MPI_Finalize(Ierror)
669 stop "Error in double torsional parameter file"
671 read (itordp,*) ntermd_1(i,j,k,iblock),
672 & ntermd_2(i,j,k,iblock)
673 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
674 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
675 read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,
676 & ntermd_1(i,j,k,iblock))
677 read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,
678 & ntermd_1(i,j,k,iblock))
679 read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,
680 & ntermd_1(i,j,k,iblock))
681 read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,
682 & ntermd_1(i,j,k,iblock))
683 C Martix of D parameters for one dimesional foureir series
684 do l=1,ntermd_1(i,j,k,iblock)
685 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
686 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
687 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
688 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
689 c write(iout,*) "whcodze" ,
690 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
692 read (itordp,*) ((v2c(l,m,i,j,k,iblock),
693 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
694 & v2s(m,l,i,j,k,iblock),
695 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
696 C Martix of D parameters for two dimesional fourier series
697 do l=1,ntermd_2(i,j,k,iblock)
699 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
700 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
701 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
702 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
711 write (iout,*) 'Constants for double torsionals'
714 do j=-ntortyp+1,ntortyp-1
715 do k=-ntortyp+1,ntortyp-1
716 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
717 & ' nsingle',ntermd_1(i,j,k,iblock),
718 & ' ndouble',ntermd_2(i,j,k,iblock)
720 write (iout,*) 'Single angles:'
721 do l=1,ntermd_1(i,j,k,iblock)
722 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
723 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
724 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
725 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
728 write (iout,*) 'Pairs of angles:'
729 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
730 do l=1,ntermd_2(i,j,k,iblock)
731 write (iout,'(i5,20f10.5)')
732 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
735 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
736 do l=1,ntermd_2(i,j,k,iblock)
737 write (iout,'(i5,20f10.5)')
738 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
739 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
748 C Read of Side-chain backbone correlation parameters
749 C Modified 11 May 2012 by Adasko
752 read (isccor,*) nsccortyp
753 read (isccor,*) (isccortyp(i),i=1,ntyp)
755 isccortyp(i)=-isccortyp(-i)
757 iscprol=isccortyp(20)
758 c write (iout,*) 'ntortyp',ntortyp
760 cc maxinter is maximum interaction sites
765 &nterm_sccor(i,j),nlor_sccor(i,j)
766 write (iout,*) nterm_sccor(i,j)
772 nterm_sccor(-i,j)=nterm_sccor(i,j)
773 nterm_sccor(-i,-j)=nterm_sccor(i,j)
774 nterm_sccor(i,-j)=nterm_sccor(i,j)
775 write (iout,*) nterm_sccor(i,j),nterm_sccor(-i,j),
776 & nterm_sccor(-i,-j),nterm_sccor(i,-j)
777 do k=1,nterm_sccor(i,j)
778 read (isccor,*) kk,v1sccor(k,l,i,j)
780 if (j.eq.iscprol) then
781 if (i.eq.isccortyp(10)) then
782 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
783 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
785 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
786 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
787 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
788 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
789 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
790 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
791 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
792 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
795 if (i.eq.isccortyp(10)) then
796 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
797 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
799 if (j.eq.isccortyp(10)) then
800 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
801 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
803 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
804 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
805 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
806 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
807 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
808 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
812 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
813 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
814 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
815 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
818 do k=1,nlor_sccor(i,j)
819 read (isccor,*) kk,vlor1sccor(k,i,j),
820 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
821 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
822 &(1+vlor3sccor(k,i,j)**2)
824 v0sccor(l,i,j)=v0ijsccor
825 v0sccor(l,-i,j)=v0ijsccor1
826 v0sccor(l,i,-j)=v0ijsccor2
827 v0sccor(l,-i,-j)=v0ijsccor3
833 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
836 write (iout,*) 'ityp',i,' jtyp',j
837 write (iout,*) 'Fourier constants'
838 do k=1,nterm_sccor(i,j)
839 write (iout,'(2(1pe15.5))')
840 & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
842 write (iout,*) 'Lorenz constants'
843 do k=1,nlor_sccor(i,j)
844 write (iout,'(3(1pe15.5))')
845 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
851 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
852 C interaction energy of the Gly, Ala, and Pro prototypes.
854 read (ifourier,*) nloctyp
857 read (ifourier,*) (b(ii,i),ii=1,13)
859 write (iout,*) 'Type',i
860 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
868 B1tilde(1,i) = b(3,i)
869 B1tilde(2,i) =-b(5,i)
870 B1tilde(1,-i) =-b(3,i)
871 B1tilde(2,-i) =b(5,i)
895 Ctilde(2,1,i)=-b(9,i)
897 Ctilde(1,1,-i)=b(7,i)
898 Ctilde(1,2,-i)=-b(9,i)
899 Ctilde(2,1,-i)=b(9,i)
900 Ctilde(2,2,-i)=b(7,i)
902 c Ctilde(1,1,i)=0.0d0
903 c Ctilde(1,2,i)=0.0d0
904 c Ctilde(2,1,i)=0.0d0
905 c Ctilde(2,2,i)=0.0d0
920 Dtilde(2,1,i)=-b(8,i)
922 Dtilde(1,1,-i)=b(6,i)
923 Dtilde(1,2,-i)=-b(8,i)
924 Dtilde(2,1,-i)=b(8,i)
925 Dtilde(2,2,-i)=b(6,i)
927 c Dtilde(1,1,i)=0.0d0
928 c Dtilde(1,2,i)=0.0d0
929 c Dtilde(2,1,i)=0.0d0
930 c Dtilde(2,2,i)=0.0d0
931 EE(1,1,i)= b(10,i)+b(11,i)
932 EE(2,2,i)=-b(10,i)+b(11,i)
933 EE(2,1,i)= b(12,i)-b(13,i)
934 EE(1,2,i)= b(12,i)+b(13,i)
935 EE(1,1,-i)= b(10,i)+b(11,i)
936 EE(2,2,-i)=-b(10,i)+b(11,i)
937 EE(2,1,-i)=-b(12,i)+b(13,i)
938 EE(1,2,-i)=-b(12,i)-b(13,i)
944 c ee(2,1,i)=ee(1,2,i)
949 write (iout,*) 'Type',i
951 c write (iout,'(f10.5)') B1(:,i)
952 write(iout,*) B1(1,i),B1(2,i)
954 c write (iout,'(f10.5)') B2(:,i)
955 write(iout,*) B2(1,i),B2(2,i)
958 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
962 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
966 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
971 C Read electrostatic-interaction parameters
974 write (iout,'(/a)') 'Electrostatic interaction constants:'
975 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
976 & 'IT','JT','APP','BPP','AEL6','AEL3'
978 read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
979 read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
980 read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
981 read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
986 app (i,j)=epp(i,j)*rri*rri
987 bpp (i,j)=-2.0D0*epp(i,j)*rri
988 ael6(i,j)=elpp6(i,j)*4.2D0**6
989 ael3(i,j)=elpp3(i,j)*4.2D0**3
991 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
992 & ael6(i,j),ael3(i,j)
997 C Read side-chain interaction parameters.
999 read (isidep,*) ipot,expon
1000 if (ipot.lt.1 .or. ipot.gt.5) then
1001 write (iout,'(2a)') 'Error while reading SC interaction',
1002 & 'potential file - unknown potential type.'
1006 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1007 & ', exponents are ',expon,2*expon
1008 goto (10,20,30,30,40) ipot
1009 C----------------------- LJ potential ---------------------------------
1010 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1012 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1013 write (iout,'(a/)') 'The epsilon array:'
1014 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1015 write (iout,'(/a)') 'One-body parameters:'
1016 write (iout,'(a,4x,a)') 'residue','sigma'
1017 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1020 C----------------------- LJK potential --------------------------------
1021 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1022 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1024 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1025 write (iout,'(a/)') 'The epsilon array:'
1026 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1027 write (iout,'(/a)') 'One-body parameters:'
1028 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1029 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1033 C---------------------- GB or BP potential -----------------------------
1035 read (isidep,*)(eps(i,j),j=i,ntyp)
1037 read (isidep,*)(sigma0(i),i=1,ntyp)
1038 read (isidep,*)(sigii(i),i=1,ntyp)
1039 read (isidep,*)(chip(i),i=1,ntyp)
1040 read (isidep,*)(alp(i),i=1,ntyp)
1042 read (isidep,*)(epslip(i,j),j=i,ntyp)
1043 C print *,"WARNING!!"
1045 C epslip(i,j)=epslip(i,j)+0.05d0
1048 C For the GB potential convert sigma'**2 into chi'
1051 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1055 write (iout,'(/a/)') 'Parameters of the BP potential:'
1056 write (iout,'(a/)') 'The epsilon array:'
1057 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1058 write (iout,'(/a)') 'One-body parameters:'
1059 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1061 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1062 & chip(i),alp(i),i=1,ntyp)
1065 C--------------------- GBV potential -----------------------------------
1066 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1067 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1068 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1070 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1071 write (iout,'(a/)') 'The epsilon array:'
1072 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1073 write (iout,'(/a)') 'One-body parameters:'
1074 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1075 & 's||/s_|_^2',' chip ',' alph '
1076 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1077 & sigii(i),chip(i),alp(i),i=1,ntyp)
1081 C-----------------------------------------------------------------------
1082 C Calculate the "working" parameters of SC interactions.
1086 epslip(i,j)=epslip(j,i)
1091 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1092 sigma(j,i)=sigma(i,j)
1093 rs0(i,j)=dwa16*sigma(i,j)
1097 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1098 & 'Working parameters of the SC interactions:',
1099 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1104 epsijlip=epslip(i,j)
1105 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1114 sigeps=dsign(1.0D0,epsij)
1116 aa_aq(i,j)=epsij*rrij*rrij
1117 bb_aq(i,j)=-sigeps*epsij*rrij
1118 aa_aq(j,i)=aa_aq(i,j)
1119 bb_aq(j,i)=bb_aq(i,j)
1120 sigeps=dsign(1.0D0,epsijlip)
1121 epsijlip=dabs(epsijlip)
1122 aa_lip(i,j)=epsijlip*rrij*rrij
1123 bb_lip(i,j)=-sigeps*epsijlip*rrij
1124 aa_lip(j,i)=aa_lip(i,j)
1125 bb_lip(j,i)=bb_lip(i,j)
1127 sigt1sq=sigma0(i)**2
1128 sigt2sq=sigma0(j)**2
1131 ratsig1=sigt2sq/sigt1sq
1132 ratsig2=1.0D0/ratsig1
1133 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1134 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1135 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1139 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1140 sigmaii(i,j)=rsum_max
1141 sigmaii(j,i)=rsum_max
1143 c sigmaii(i,j)=r0(i,j)
1144 c sigmaii(j,i)=r0(i,j)
1146 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1147 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1148 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1149 augm(i,j)=epsij*r_augm**(2*expon)
1150 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1157 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1158 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1159 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1164 C Define the SC-p interaction constants
1168 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1170 c aad(i,1)=0.3D0*4.0D0**12
1171 C Following line for constants currently implemented
1172 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1173 aad(i,1)=1.5D0*4.0D0**12
1174 c aad(i,1)=0.17D0*5.6D0**12
1176 C "Soft" SC-p repulsion
1178 C Following line for constants currently implemented
1179 c aad(i,1)=0.3D0*4.0D0**6
1180 C "Hard" SC-p repulsion
1181 bad(i,1)=3.0D0*4.0D0**6
1182 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1191 C 8/9/01 Read the SC-p interaction constants from file
1194 read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1197 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1198 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1199 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1200 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1204 write (iout,*) "Parameters of SC-p interactions:"
1206 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1207 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1212 C Define the constants of the disulfide bridge
1216 c Old arbitrary potential - commented out.
1221 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1222 c energy surface of diethyl disulfide.
1223 c A. Liwo and U. Kozlowska, 11/24/03
1234 write (iout,'(/a)') "Disulfide bridge parameters:"
1235 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1236 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1237 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1238 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,