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)
545 itortyp(i)=-itortyp(-i)
547 c write (iout,*) 'ntortyp',ntortyp
549 do j=-ntortyp+1,ntortyp-1
550 read (itorp,*) nterm(i,j,iblock),
552 nterm(-i,-j,iblock)=nterm(i,j,iblock)
553 nlor(-i,-j,iblock)=nlor(i,j,iblock)
556 do k=1,nterm(i,j,iblock)
557 read (itorp,*) kk,v1(k,i,j,iblock),v2(k,i,j,iblock)
558 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
559 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
560 v0ij=v0ij+si*v1(k,i,j,iblock)
563 do k=1,nlor(i,j,iblock)
564 read (itorp,*) kk,vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
565 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
568 v0(-i,-j,iblock)=v0ij
574 write (iout,'(/a/)') 'Torsional constants:'
577 write (iout,*) 'ityp',i,' jtyp',j
578 write (iout,*) 'Fourier constants'
579 do k=1,nterm(i,j,iblock)
580 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
583 write (iout,*) 'Lorenz constants'
584 do k=1,nlor(i,j,iblock)
585 write (iout,'(3(1pe15.5))')
586 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
592 C 6/23/01 Read parameters for double torsionals
596 do j=-ntortyp+1,ntortyp-1
597 do k=-ntortyp+1,ntortyp-1
598 read (itordp,'(3a1)') t1,t2,t3
599 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
600 & .or. t3.ne.toronelet(k)) then
601 write (iout,*) "Error in double torsional parameter file",
603 stop "Error in double torsional parameter file"
605 read (itordp,*) ntermd_1(i,j,k,iblock),
606 & ntermd_2(i,j,k,iblock)
607 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
608 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
609 read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,
610 & ntermd_1(i,j,k,iblock))
611 read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,
612 & ntermd_1(i,j,k,iblock))
613 read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,
614 & ntermd_1(i,j,k,iblock))
615 read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,
616 & ntermd_1(i,j,k,iblock))
617 C Martix of D parameters for one dimesional foureir series
618 do l=1,ntermd_1(i,j,k,iblock)
619 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
620 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
621 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
622 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
623 c write(iout,*) "whcodze" ,
624 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
626 read (itordp,*) ((v2c(l,m,i,j,k,iblock),
627 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
628 & v2s(m,l,i,j,k,iblock),
629 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
630 C Martix of D parameters for two dimesional fourier series
631 do l=1,ntermd_2(i,j,k,iblock)
633 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
634 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
635 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
636 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
645 write (iout,*) 'Constants for double torsionals'
648 do j=-ntortyp+1,ntortyp-1
649 do k=-ntortyp+1,ntortyp-1
650 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
651 & ' nsingle',ntermd_1(i,j,k,iblock),
652 & ' ndouble',ntermd_2(i,j,k,iblock)
654 write (iout,*) 'Single angles:'
655 do l=1,ntermd_1(i,j,k,iblock)
656 write (iout,'(i5,2f10.5,5x,2f10.5)') l,
657 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
658 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock)
661 write (iout,*) 'Pairs of angles:'
662 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
663 do l=1,ntermd_2(i,j,k,iblock)
664 write (iout,'(i5,20f10.5)')
665 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
668 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
669 do l=1,ntermd_2(i,j,k,iblock)
670 write (iout,'(i5,20f10.5)')
671 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
680 C Read of Side-chain backbone correlation parameters
681 C Modified 11 May 2012 by Adasko
684 read (isccor,*) nsccortyp
685 read (isccor,*) (isccortyp(i),i=1,ntyp)
687 isccortyp(i)=-isccortyp(-i)
689 iscprol=isccortyp(20)
690 c write (iout,*) 'ntortyp',ntortyp
692 cc maxinter is maximum interaction sites
696 read (isccor,*) nterm_sccor(i,j),nlor_sccor(i,j)
702 nterm_sccor(-i,j)=nterm_sccor(i,j)
703 nterm_sccor(-i,-j)=nterm_sccor(i,j)
704 nterm_sccor(i,-j)=nterm_sccor(i,j)
705 do k=1,nterm_sccor(i,j)
706 read (isccor,*) kk,v1sccor(k,l,i,j)
708 if (j.eq.iscprol) then
709 if (i.eq.isccortyp(10)) then
710 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
711 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
713 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
714 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
715 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
716 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
717 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
718 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
719 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
720 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
723 if (i.eq.isccortyp(10)) then
724 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
725 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
727 if (j.eq.isccortyp(10)) then
728 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
729 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
731 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
732 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
733 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
734 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
735 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
736 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
740 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
741 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
742 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
743 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
746 do k=1,nlor_sccor(i,j)
747 read (isccor,*) kk,vlor1sccor(k,i,j),
748 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
749 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
750 &(1+vlor3sccor(k,i,j)**2)
752 v0sccor(l,i,j)=v0ijsccor
753 v0sccor(l,-i,j)=v0ijsccor1
754 v0sccor(l,i,-j)=v0ijsccor2
755 v0sccor(l,-i,-j)=v0ijsccor3
762 write (iout,'(/a/)') 'Torsional constants:'
765 write (iout,*) 'ityp',i,' jtyp',j
766 write (iout,*) 'Fourier constants'
767 do k=1,nterm_sccor(i,j)
768 write (iout,'(2(1pe15.5))')v1sccor(k,l,i,j),v2sccor(k,l,i,j)
770 write (iout,*) 'Lorenz constants'
771 do k=1,nlor_sccor(i,j)
772 write (iout,'(3(1pe15.5))')
773 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
780 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
781 C interaction energy of the Gly, Ala, and Pro prototypes.
783 read (ifourier,*) nloctyp
786 read (ifourier,*) (b(ii),ii=1,13)
788 write (iout,*) 'Type',i
789 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
831 c Ctilde(1,1,i)=0.0d0
832 c Ctilde(1,2,i)=0.0d0
833 c Ctilde(2,1,i)=0.0d0
834 c Ctilde(2,2,i)=0.0d0
856 c Dtilde(1,1,i)=0.0d0
857 c Dtilde(1,2,i)=0.0d0
858 c Dtilde(2,1,i)=0.0d0
859 c Dtilde(2,2,i)=0.0d0
860 EE(1,1,i)= b(10)+b(11)
861 EE(2,2,i)=-b(10)+b(11)
862 EE(2,1,i)= b(12)-b(13)
863 EE(1,2,i)= b(12)+b(13)
864 EE(1,1,-i)= b(10)+b(11)
865 EE(2,2,-i)=-b(10)+b(11)
866 EE(2,1,-i)=-b(12)+b(13)
867 EE(1,2,-i)=-b(12)-b(13)
873 c ee(2,1,i)=ee(1,2,i)
877 write (iout,*) 'Type',i
879 c write (iout,'(f10.5)') B1(:,i)
880 write(iout,*) B1(1,i),B1(2,i)
882 c write (iout,'(f10.5)') B2(:,i)
883 write(iout,*) B2(1,i),B2(2,i)
886 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
890 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
894 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
899 C Read electrostatic-interaction parameters
902 write (iout,'(/a)') 'Electrostatic interaction constants:'
903 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
904 & 'IT','JT','APP','BPP','AEL6','AEL3'
906 read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
907 read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
908 read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
909 read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
914 app (i,j)=epp(i,j)*rri*rri
915 bpp (i,j)=-2.0D0*epp(i,j)*rri
916 ael6(i,j)=elpp6(i,j)*4.2D0**6
917 ael3(i,j)=elpp3(i,j)*4.2D0**3
918 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
919 & ael6(i,j),ael3(i,j)
923 C Read side-chain interaction parameters.
925 read (isidep,*) ipot,expon
926 if (ipot.lt.1 .or. ipot.gt.5) then
927 write (iout,'(2a)') 'Error while reading SC interaction',
928 & 'potential file - unknown potential type.'
932 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
933 & ', exponents are ',expon,2*expon
934 goto (10,20,30,30,40) ipot
935 C----------------------- LJ potential ---------------------------------
936 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
938 write (iout,'(/a/)') 'Parameters of the LJ potential:'
939 write (iout,'(a/)') 'The epsilon array:'
940 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
941 write (iout,'(/a)') 'One-body parameters:'
942 write (iout,'(a,4x,a)') 'residue','sigma'
943 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
946 C----------------------- LJK potential --------------------------------
947 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
948 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
950 write (iout,'(/a/)') 'Parameters of the LJK potential:'
951 write (iout,'(a/)') 'The epsilon array:'
952 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
953 write (iout,'(/a)') 'One-body parameters:'
954 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
955 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
959 C---------------------- GB or BP potential -----------------------------
960 30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
961 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
963 C For the GB potential convert sigma'**2 into chi'
966 chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
970 write (iout,'(/a/)') 'Parameters of the BP potential:'
971 write (iout,'(a/)') 'The epsilon array:'
972 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
973 write (iout,'(/a)') 'One-body parameters:'
974 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
976 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
977 & chip(i),alp(i),i=1,ntyp)
980 C--------------------- GBV potential -----------------------------------
981 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
982 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
983 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
985 write (iout,'(/a/)') 'Parameters of the GBV potential:'
986 write (iout,'(a/)') 'The epsilon array:'
987 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
988 write (iout,'(/a)') 'One-body parameters:'
989 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
990 & 's||/s_|_^2',' chip ',' alph '
991 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
992 & sigii(i),chip(i),alp(i),i=1,ntyp)
996 C-----------------------------------------------------------------------
997 C Calculate the "working" parameters of SC interactions.
1005 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1006 sigma(j,i)=sigma(i,j)
1007 rs0(i,j)=dwa16*sigma(i,j)
1011 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1012 & 'Working parameters of the SC interactions:',
1013 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1018 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1027 sigeps=dsign(1.0D0,epsij)
1029 aa(i,j)=epsij*rrij*rrij
1030 bb(i,j)=-sigeps*epsij*rrij
1034 sigt1sq=sigma0(i)**2
1035 sigt2sq=sigma0(j)**2
1038 ratsig1=sigt2sq/sigt1sq
1039 ratsig2=1.0D0/ratsig1
1040 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1041 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1042 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1046 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1047 sigmaii(i,j)=rsum_max
1048 sigmaii(j,i)=rsum_max
1050 c sigmaii(i,j)=r0(i,j)
1051 c sigmaii(j,i)=r0(i,j)
1053 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1054 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1055 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1056 augm(i,j)=epsij*r_augm**(2*expon)
1057 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1064 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1065 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1066 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1071 C Define the SC-p interaction constants
1075 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1077 c aad(i,1)=0.3D0*4.0D0**12
1078 C Following line for constants currently implemented
1079 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1080 aad(i,1)=1.5D0*4.0D0**12
1081 c aad(i,1)=0.17D0*5.6D0**12
1083 C "Soft" SC-p repulsion
1085 C Following line for constants currently implemented
1086 c aad(i,1)=0.3D0*4.0D0**6
1087 C "Hard" SC-p repulsion
1088 bad(i,1)=3.0D0*4.0D0**6
1089 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1098 C 8/9/01 Read the SC-p interaction constants from file
1101 read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1104 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1105 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1106 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1107 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1111 write (iout,*) "Parameters of SC-p interactions:"
1113 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1114 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1119 C Define the constants of the disulfide bridge
1123 c Old arbitrary potential - commented out.
1128 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1129 c energy surface of diethyl disulfide.
1130 c A. Liwo and U. Kozlowska, 11/24/03
1141 write (iout,'(/a)') "Disulfide bridge parameters:"
1142 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1143 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1144 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1145 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,