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 write (iout,*) "tu dochodze"
299 read (ithep,*) nthetyp,ntheterm,ntheterm2,
300 & ntheterm3,nsingle,ndouble
301 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
302 read (ithep,*) (ithetyp(i),i=1,ntyp1)
304 ithetyp(i)=-ithetyp(-i)
306 write (iout,*) "tu dochodze"
308 do i=-maxthetyp,maxthetyp
309 do j=-maxthetyp,maxthetyp
310 do k=-maxthetyp,maxthetyp
311 aa0thet(i,j,k,iblock)=0.0d0
313 aathet(l,i,j,k,iblock)=0.0d0
317 bbthet(m,l,i,j,k,iblock)=0.0d0
318 ccthet(m,l,i,j,k,iblock)=0.0d0
319 ddthet(m,l,i,j,k,iblock)=0.0d0
320 eethet(m,l,i,j,k,iblock)=0.0d0
326 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
327 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
337 do j=-nthetyp,nthetyp
338 do k=-nthetyp,nthetyp
339 read (ithep,'(6a)') res1
340 read (ithep,*) aa0thet(i,j,k,iblock)
341 read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
343 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
344 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
345 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
346 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle)
349 & (((ffthet(llll,lll,ll,i,j,k,iblock),
350 & ffthet(lll,llll,ll,i,j,k,iblock),
351 & ggthet(llll,lll,ll,i,j,k,iblock)
352 & ,ggthet(lll,llll,ll,i,j,k,iblock),
353 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
358 C For dummy ends assign glycine-type coefficients of theta-only terms; the
359 C coefficients of theta-and-gamma-dependent terms are zero.
364 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
365 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
367 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
368 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
371 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
373 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
376 C Substitution for D aminoacids from symmetry.
379 do j=-nthetyp,nthetyp
380 do k=-nthetyp,nthetyp
381 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
383 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
387 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
388 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
389 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
390 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
396 ffthet(llll,lll,ll,i,j,k,iblock)=
397 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
398 ffthet(lll,llll,ll,i,j,k,iblock)=
399 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
400 ggthet(llll,lll,ll,i,j,k,iblock)=
401 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
402 ggthet(lll,llll,ll,i,j,k,iblock)=
403 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
413 C Control printout of the coefficients of virtual-bond-angle potentials
416 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
420 write (iout,'(//4a)')
421 & 'Type ',onelett(i),onelett(j),onelett(k)
422 write (iout,'(//a,10x,a)') " l","a[l]"
423 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
424 write (iout,'(i2,1pe15.5)')
425 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
427 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
428 & "b",l,"c",l,"d",l,"e",l
430 write (iout,'(i2,4(1pe15.5))') m,
431 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
432 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
436 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
437 & "f+",l,"f-",l,"g+",l,"g-",l
440 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
441 & ffthet(n,m,l,i,j,k,iblock),
442 & ffthet(m,n,l,i,j,k,iblock),
443 & ggthet(n,m,l,i,j,k,iblock),
444 & ggthet(m,n,l,i,j,k,iblock)
457 C Read the parameters of the probability distribution/energy expression
458 C of the side chains.
461 cc write (iout,*) "tu dochodze",i
462 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
466 dsc_inv(i)=1.0D0/dsc(i)
477 read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
478 censc(1,1,-i)=censc(1,1,i)
479 censc(2,1,-i)=censc(2,1,i)
480 censc(3,1,-i)=-censc(3,1,i)
482 read (irotam,*) bsc(j,i)
483 read (irotam,*) (censc(k,j,i),k=1,3),
484 & ((blower(k,l,j),l=1,k),k=1,3)
485 censc(1,j,-i)=censc(1,j,i)
486 censc(2,j,-i)=censc(2,j,i)
487 censc(3,j,-i)=-censc(3,j,i)
488 C BSC is amplitude of Gaussian
495 akl=akl+blower(k,m,j)*blower(l,m,j)
499 if (((k.eq.3).and.(l.ne.3))
500 & .or.((l.eq.3).and.(k.ne.3))) then
501 gaussc(k,l,j,-i)=-akl
502 gaussc(l,k,j,-i)=-akl
514 write (iout,'(/a)') 'Parameters of side-chain local geometry'
518 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
519 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
520 c write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
521 c write (iout,'(a,f10.4,4(16x,f10.4))')
522 c & 'Center ',(bsc(j,i),j=1,nlobi)
523 c write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
524 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
525 & 'log h',(bsc(j,i),j=1,nlobi)
526 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
527 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
534 c blower(k,l,j)=gaussc(ind,j,i)
539 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
540 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
547 C Read scrot parameters for potentials determined from all-atom AM1 calculations
548 C added by Urszula Kozlowska 07/11/2007
556 read(irotam,*) sc_parmin(j,i)
564 C Read torsional parameters in old format
566 read (itorp,*) ntortyp,nterm_old
567 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
568 read (itorp,*) (itortyp(i),i=1,ntyp)
573 read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
579 write (iout,'(/a/)') 'Torsional constants:'
582 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
583 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
591 C Read torsional parameters
593 read (itorp,*) ntortyp
594 read (itorp,*) (itortyp(i),i=1,ntyp)
595 write (iout,*) 'ntortyp',ntortyp
598 itortyp(i)=-itortyp(-i)
600 c write (iout,*) 'ntortyp',ntortyp
602 do j=-ntortyp+1,ntortyp-1
603 read (itorp,*) nterm(i,j,iblock),
605 nterm(-i,-j,iblock)=nterm(i,j,iblock)
606 nlor(-i,-j,iblock)=nlor(i,j,iblock)
609 do k=1,nterm(i,j,iblock)
610 read (itorp,*) kk,v1(k,i,j,iblock),
612 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
613 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
614 v0ij=v0ij+si*v1(k,i,j,iblock)
617 do k=1,nlor(i,j,iblock)
618 read (itorp,*) kk,vlor1(k,i,j),
619 & vlor2(k,i,j),vlor3(k,i,j)
620 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
623 v0(-i,-j,iblock)=v0ij
629 write (iout,'(/a/)') 'Torsional constants:'
632 write (iout,*) 'ityp',i,' jtyp',j
633 write (iout,*) 'Fourier constants'
634 do k=1,nterm(i,j,iblock)
635 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
638 write (iout,*) 'Lorenz constants'
639 do k=1,nlor(i,j,iblock)
640 write (iout,'(3(1pe15.5))')
641 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
647 C 6/23/01 Read parameters for double torsionals
651 do j=-ntortyp+1,ntortyp-1
652 do k=-ntortyp+1,ntortyp-1
653 read (itordp,'(3a1)') t1,t2,t3
654 c write (iout,*) "OK onelett",
657 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
658 & .or. t3.ne.toronelet(k)) then
659 write (iout,*) "Error in double torsional parameter file",
662 call MPI_Finalize(Ierror)
664 stop "Error in double torsional parameter file"
666 read (itordp,*) ntermd_1(i,j,k,iblock),
667 & ntermd_2(i,j,k,iblock)
668 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
669 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
670 read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,
671 & ntermd_1(i,j,k,iblock))
672 read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,
673 & ntermd_1(i,j,k,iblock))
674 read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,
675 & ntermd_1(i,j,k,iblock))
676 read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,
677 & ntermd_1(i,j,k,iblock))
678 C Martix of D parameters for one dimesional foureir series
679 do l=1,ntermd_1(i,j,k,iblock)
680 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
681 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
682 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
683 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
684 c write(iout,*) "whcodze" ,
685 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
687 read (itordp,*) ((v2c(l,m,i,j,k,iblock),
688 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
689 & v2s(m,l,i,j,k,iblock),
690 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
691 C Martix of D parameters for two dimesional fourier series
692 do l=1,ntermd_2(i,j,k,iblock)
694 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
695 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
696 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
697 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
706 write (iout,*) 'Constants for double torsionals'
709 do j=-ntortyp+1,ntortyp-1
710 do k=-ntortyp+1,ntortyp-1
711 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
712 & ' nsingle',ntermd_1(i,j,k,iblock),
713 & ' ndouble',ntermd_2(i,j,k,iblock)
715 write (iout,*) 'Single angles:'
716 do l=1,ntermd_1(i,j,k,iblock)
717 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
718 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
719 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
720 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
723 write (iout,*) 'Pairs of angles:'
724 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
725 do l=1,ntermd_2(i,j,k,iblock)
726 write (iout,'(i5,20f10.5)')
727 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
730 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
731 do l=1,ntermd_2(i,j,k,iblock)
732 write (iout,'(i5,20f10.5)')
733 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
734 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
743 C Read of Side-chain backbone correlation parameters
744 C Modified 11 May 2012 by Adasko
747 read (isccor,*) nsccortyp
748 read (isccor,*) (isccortyp(i),i=1,ntyp)
750 isccortyp(i)=-isccortyp(-i)
752 iscprol=isccortyp(20)
753 c write (iout,*) 'ntortyp',ntortyp
755 cc maxinter is maximum interaction sites
760 &nterm_sccor(i,j),nlor_sccor(i,j)
761 write (iout,*) nterm_sccor(i,j)
767 nterm_sccor(-i,j)=nterm_sccor(i,j)
768 nterm_sccor(-i,-j)=nterm_sccor(i,j)
769 nterm_sccor(i,-j)=nterm_sccor(i,j)
770 write (iout,*) nterm_sccor(i,j),nterm_sccor(-i,j),
771 & nterm_sccor(-i,-j),nterm_sccor(i,-j)
772 do k=1,nterm_sccor(i,j)
773 read (isccor,*) kk,v1sccor(k,l,i,j)
775 if (j.eq.iscprol) then
776 if (i.eq.isccortyp(10)) then
777 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
778 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
780 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
781 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
782 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
783 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
784 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
785 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
786 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
787 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
790 if (i.eq.isccortyp(10)) then
791 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
792 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
794 if (j.eq.isccortyp(10)) then
795 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
796 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)
802 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
803 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
807 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
808 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
809 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
810 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
813 do k=1,nlor_sccor(i,j)
814 read (isccor,*) kk,vlor1sccor(k,i,j),
815 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
816 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
817 &(1+vlor3sccor(k,i,j)**2)
819 v0sccor(l,i,j)=v0ijsccor
820 v0sccor(l,-i,j)=v0ijsccor1
821 v0sccor(l,i,-j)=v0ijsccor2
822 v0sccor(l,-i,-j)=v0ijsccor3
828 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
831 write (iout,*) 'ityp',i,' jtyp',j
832 write (iout,*) 'Fourier constants'
833 do k=1,nterm_sccor(i,j)
834 write (iout,'(2(1pe15.5))')
835 & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
837 write (iout,*) 'Lorenz constants'
838 do k=1,nlor_sccor(i,j)
839 write (iout,'(3(1pe15.5))')
840 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
846 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
847 C interaction energy of the Gly, Ala, and Pro prototypes.
849 read (ifourier,*) nloctyp
852 read (ifourier,*) (b(ii,i),ii=1,13)
854 write (iout,*) 'Type',i
855 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
863 B1tilde(1,i) = b(3,i)
864 B1tilde(2,i) =-b(5,i)
865 B1tilde(1,-i) =-b(3,i)
866 B1tilde(2,-i) =b(5,i)
890 Ctilde(2,1,i)=-b(9,i)
892 Ctilde(1,1,-i)=b(7,i)
893 Ctilde(1,2,-i)=-b(9,i)
894 Ctilde(2,1,-i)=b(9,i)
895 Ctilde(2,2,-i)=b(7,i)
897 c Ctilde(1,1,i)=0.0d0
898 c Ctilde(1,2,i)=0.0d0
899 c Ctilde(2,1,i)=0.0d0
900 c Ctilde(2,2,i)=0.0d0
915 Dtilde(2,1,i)=-b(8,i)
917 Dtilde(1,1,-i)=b(6,i)
918 Dtilde(1,2,-i)=-b(8,i)
919 Dtilde(2,1,-i)=b(8,i)
920 Dtilde(2,2,-i)=b(6,i)
922 c Dtilde(1,1,i)=0.0d0
923 c Dtilde(1,2,i)=0.0d0
924 c Dtilde(2,1,i)=0.0d0
925 c Dtilde(2,2,i)=0.0d0
926 EE(1,1,i)= b(10,i)+b(11,i)
927 EE(2,2,i)=-b(10,i)+b(11,i)
928 EE(2,1,i)= b(12,i)-b(13,i)
929 EE(1,2,i)= b(12,i)+b(13,i)
930 EE(1,1,-i)= b(10,i)+b(11,i)
931 EE(2,2,-i)=-b(10,i)+b(11,i)
932 EE(2,1,-i)=-b(12,i)+b(13,i)
933 EE(1,2,-i)=-b(12,i)-b(13,i)
939 c ee(2,1,i)=ee(1,2,i)
944 write (iout,*) 'Type',i
946 c write (iout,'(f10.5)') B1(:,i)
947 write(iout,*) B1(1,i),B1(2,i)
949 c write (iout,'(f10.5)') B2(:,i)
950 write(iout,*) B2(1,i),B2(2,i)
953 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
957 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
961 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
966 C Read electrostatic-interaction parameters
969 write (iout,'(/a)') 'Electrostatic interaction constants:'
970 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
971 & 'IT','JT','APP','BPP','AEL6','AEL3'
973 read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
974 read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
975 read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
976 read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
981 app (i,j)=epp(i,j)*rri*rri
982 bpp (i,j)=-2.0D0*epp(i,j)*rri
983 ael6(i,j)=elpp6(i,j)*4.2D0**6
984 ael3(i,j)=elpp3(i,j)*4.2D0**3
985 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
986 & ael6(i,j),ael3(i,j)
990 C Read side-chain interaction parameters.
992 read (isidep,*) ipot,expon
993 if (ipot.lt.1 .or. ipot.gt.5) then
994 write (iout,'(2a)') 'Error while reading SC interaction',
995 & 'potential file - unknown potential type.'
999 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1000 & ', exponents are ',expon,2*expon
1001 goto (10,20,30,30,40) ipot
1002 C----------------------- LJ potential ---------------------------------
1003 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1005 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1006 write (iout,'(a/)') 'The epsilon array:'
1007 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1008 write (iout,'(/a)') 'One-body parameters:'
1009 write (iout,'(a,4x,a)') 'residue','sigma'
1010 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1013 C----------------------- LJK potential --------------------------------
1014 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1015 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1017 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1018 write (iout,'(a/)') 'The epsilon array:'
1019 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1020 write (iout,'(/a)') 'One-body parameters:'
1021 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1022 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1026 C---------------------- GB or BP potential -----------------------------
1028 read (isidep,*)(eps(i,j),j=i,ntyp)
1030 read (isidep,*)(sigma0(i),i=1,ntyp)
1031 read (isidep,*)(sigii(i),i=1,ntyp)
1032 read (isidep,*)(chip(i),i=1,ntyp)
1033 read (isidep,*)(alp(i),i=1,ntyp)
1034 C For the GB potential convert sigma'**2 into chi'
1037 chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
1041 write (iout,'(/a/)') 'Parameters of the BP potential:'
1042 write (iout,'(a/)') 'The epsilon array:'
1043 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1044 write (iout,'(/a)') 'One-body parameters:'
1045 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1047 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1048 & chip(i),alp(i),i=1,ntyp)
1051 C--------------------- GBV potential -----------------------------------
1052 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1053 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1054 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1056 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1057 write (iout,'(a/)') 'The epsilon array:'
1058 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1059 write (iout,'(/a)') 'One-body parameters:'
1060 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1061 & 's||/s_|_^2',' chip ',' alph '
1062 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1063 & sigii(i),chip(i),alp(i),i=1,ntyp)
1067 C-----------------------------------------------------------------------
1068 C Calculate the "working" parameters of SC interactions.
1076 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1077 sigma(j,i)=sigma(i,j)
1078 rs0(i,j)=dwa16*sigma(i,j)
1082 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1083 & 'Working parameters of the SC interactions:',
1084 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1089 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1098 sigeps=dsign(1.0D0,epsij)
1100 aa(i,j)=epsij*rrij*rrij
1101 bb(i,j)=-sigeps*epsij*rrij
1105 sigt1sq=sigma0(i)**2
1106 sigt2sq=sigma0(j)**2
1109 ratsig1=sigt2sq/sigt1sq
1110 ratsig2=1.0D0/ratsig1
1111 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1112 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1113 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1117 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1118 sigmaii(i,j)=rsum_max
1119 sigmaii(j,i)=rsum_max
1121 c sigmaii(i,j)=r0(i,j)
1122 c sigmaii(j,i)=r0(i,j)
1124 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1125 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1126 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1127 augm(i,j)=epsij*r_augm**(2*expon)
1128 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1135 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1136 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1137 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1142 C Define the SC-p interaction constants
1146 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1148 c aad(i,1)=0.3D0*4.0D0**12
1149 C Following line for constants currently implemented
1150 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1151 aad(i,1)=1.5D0*4.0D0**12
1152 c aad(i,1)=0.17D0*5.6D0**12
1154 C "Soft" SC-p repulsion
1156 C Following line for constants currently implemented
1157 c aad(i,1)=0.3D0*4.0D0**6
1158 C "Hard" SC-p repulsion
1159 bad(i,1)=3.0D0*4.0D0**6
1160 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1169 C 8/9/01 Read the SC-p interaction constants from file
1172 read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1175 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1176 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1177 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1178 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1182 write (iout,*) "Parameters of SC-p interactions:"
1184 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1185 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1190 C Define the constants of the disulfide bridge
1194 c Old arbitrary potential - commented out.
1199 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1200 c energy surface of diethyl disulfide.
1201 c A. Liwo and U. Kozlowska, 11/24/03
1212 write (iout,'(/a)') "Disulfide bridge parameters:"
1213 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1214 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1215 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1216 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,