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
58 call reada(controlcard,"DTRISS",dtriss,1.0D0)
59 call reada(controlcard,"ATRISS",atriss,0.3D0)
60 call reada(controlcard,"BTRISS",btriss,0.02D0)
61 call reada(controlcard,"CTRISS",ctriss,1.0D0)
62 dyn_ss=(index(controlcard,'DYN_SS').gt.0)
64 C dyn_ss_mask(i)=.false.
68 dyn_ssbond_ij(i,j)=1.0d300
71 call reada(controlcard,"HT",Ht,0.0D0)
73 ss_depth=ebr/wsc-0.25*eps(1,1)
74 Ht=Ht/wsc-0.25*eps(1,1)
82 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
85 if (iparm.eq.myparm .or. .not.separate_parset) then
88 c Setup weights for UNRES
110 call card_concat(controlcard,.false.)
112 c Return if not own parameters
114 if (iparm.ne.myparm .and. separate_parset) return
116 call reads(controlcard,"BONDPAR",bondname_t,bondname)
117 open (ibond,file=bondname_t,status='old')
119 call reads(controlcard,"THETPAR",thetname_t,thetname)
120 open (ithep,file=thetname_t,status='old')
122 call reads(controlcard,"ROTPAR",rotname_t,rotname)
123 open (irotam,file=rotname_t,status='old')
125 call reads(controlcard,"TORPAR",torname_t,torname)
126 open (itorp,file=torname_t,status='old')
128 call reads(controlcard,"TORDPAR",tordname_t,tordname)
129 open (itordp,file=tordname_t,status='old')
131 call reads(controlcard,"SCCORPAR",sccorname_t,sccorname)
132 open (isccor,file=sccorname_t,status='old')
134 call reads(controlcard,"FOURIER",fouriername_t,fouriername)
135 open (ifourier,file=fouriername_t,status='old')
137 call reads(controlcard,"ELEPAR",elename_t,elename)
138 open (ielep,file=elename_t,status='old')
140 call reads(controlcard,"SIDEPAR",sidename_t,sidename)
141 open (isidep,file=sidename_t,status='old')
143 call reads(controlcard,"SCPPAR",scpname_t,scpname)
144 open (iscpp,file=scpname_t,status='old')
146 write (iout,*) "Parameter set:",iparm
147 write (iout,*) "Energy-term weights:"
149 write (iout,'(a16,f10.5)') wname(i),ww(i)
151 write (iout,*) "Sidechain potential file : ",
152 & sidename_t(:ilen(sidename_t))
154 write (iout,*) "SCp potential file : ",
155 & scpname_t(:ilen(scpname_t))
157 write (iout,*) "Electrostatic potential file : ",
158 & elename_t(:ilen(elename_t))
159 write (iout,*) "Cumulant coefficient file : ",
160 & fouriername_t(:ilen(fouriername_t))
161 write (iout,*) "Torsional parameter file : ",
162 & torname_t(:ilen(torname_t))
163 write (iout,*) "Double torsional parameter file : ",
164 & tordname_t(:ilen(tordname_t))
165 write (iout,*) "Backbone-rotamer parameter file : ",
166 & sccorname(:ilen(sccorname))
167 write (iout,*) "Bond & inertia constant file : ",
168 & bondname_t(:ilen(bondname_t))
169 write (iout,*) "Bending parameter file : ",
170 & thetname_t(:ilen(thetname_t))
171 write (iout,*) "Rotamer parameter file : ",
172 & rotname_t(:ilen(rotname_t))
175 c Read the virtual-bond parameters, masses, and moments of inertia
176 c and Stokes' radii of the peptide group and side chains
179 read (ibond,*) vbldp0,akp
182 read (ibond,*) vbldsc0(1,i),aksc(1,i)
183 dsc(i) = vbldsc0(1,i)
187 dsc_inv(i)=1.0D0/dsc(i)
191 read (ibond,*) ijunk,vbldp0,akp,rjunk
193 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
195 dsc(i) = vbldsc0(1,i)
199 dsc_inv(i)=1.0D0/dsc(i)
204 write(iout,'(/a/)')"Force constants virtual bonds:"
205 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
207 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
209 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
210 & vbldsc0(1,i),aksc(1,i),abond0(1,i)
212 write (iout,'(13x,3f10.5)')
213 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
219 C Read the parameters of the probability distribution/energy expression
220 C of the virtual-bond valence angles theta
223 read (ithep,*) a0thet(i),(athet(j,i,1,1),j=1,2),
224 & (bthet(j,i,1,1),j=1,2)
225 read (ithep,*) (polthet(j,i),j=0,3)
226 read (ithep,*) (gthet(j,i),j=1,3)
227 read (ithep,*) theta0(i),sig0(i),sigc0(i)
231 athet(1,i,1,-1)=athet(1,i,1,1)
232 athet(2,i,1,-1)=athet(2,i,1,1)
233 bthet(1,i,1,-1)=-bthet(1,i,1,1)
234 bthet(2,i,1,-1)=-bthet(2,i,1,1)
235 athet(1,i,-1,1)=-athet(1,i,1,1)
236 athet(2,i,-1,1)=-athet(2,i,1,1)
237 bthet(1,i,-1,1)=bthet(1,i,1,1)
238 bthet(2,i,-1,1)=bthet(2,i,1,1)
242 athet(1,i,-1,-1)=athet(1,-i,1,1)
243 athet(2,i,-1,-1)=-athet(2,-i,1,1)
244 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
245 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
246 athet(1,i,-1,1)=athet(1,-i,1,1)
247 athet(2,i,-1,1)=-athet(2,-i,1,1)
248 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
249 bthet(2,i,-1,1)=bthet(2,-i,1,1)
250 athet(1,i,1,-1)=-athet(1,-i,1,1)
251 athet(2,i,1,-1)=athet(2,-i,1,1)
252 bthet(1,i,1,-1)=bthet(1,-i,1,1)
253 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
258 polthet(j,i)=polthet(j,-i)
261 gthet(j,i)=gthet(j,-i)
267 c & 'Parameters of the virtual-bond valence angles:'
268 c write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
269 c & ' ATHETA0 ',' A1 ',' A2 ',
272 c write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
273 c & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
275 c write (iout,'(/a/9x,5a/79(1h-))')
276 c & 'Parameters of the expression for sigma(theta_c):',
277 c & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
278 c & ' ALPH3 ',' SIGMA0C '
280 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
281 c & (polthet(j,i),j=0,3),sigc0(i)
283 c write (iout,'(/a/9x,5a/79(1h-))')
284 c & 'Parameters of the second gaussian:',
285 c & ' THETA0 ',' SIGMA0 ',' G1 ',
288 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
289 c & sig0(i),(gthet(j,i),j=1,3)
292 & 'Parameters of the virtual-bond valence angles:'
293 write (iout,'(/a/9x,5a/79(1h-))')
294 & 'Coefficients of expansion',
295 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
296 & ' b1*10^1 ',' b2*10^1 '
298 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
299 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
300 & (10*bthet(j,i,1,1),j=1,2)
302 write (iout,'(/a/9x,5a/79(1h-))')
303 & 'Parameters of the expression for sigma(theta_c):',
304 & ' alpha0 ',' alph1 ',' alph2 ',
305 & ' alhp3 ',' sigma0c '
307 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
308 & (polthet(j,i),j=0,3),sigc0(i)
310 write (iout,'(/a/9x,5a/79(1h-))')
311 & 'Parameters of the second gaussian:',
312 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
315 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
316 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
321 C Read the parameters of Utheta determined from ab initio surfaces
322 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
324 c write (iout,*) "tu dochodze"
325 read (ithep,*) nthetyp,ntheterm,ntheterm2,
326 & ntheterm3,nsingle,ndouble
327 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
328 read (ithep,*) (ithetyp(i),i=1,ntyp1)
330 ithetyp(i)=-ithetyp(-i)
332 c write (iout,*) "tu dochodze"
334 do i=-maxthetyp,maxthetyp
335 do j=-maxthetyp,maxthetyp
336 do k=-maxthetyp,maxthetyp
337 aa0thet(i,j,k,iblock)=0.0d0
339 aathet(l,i,j,k,iblock)=0.0d0
343 bbthet(m,l,i,j,k,iblock)=0.0d0
344 ccthet(m,l,i,j,k,iblock)=0.0d0
345 ddthet(m,l,i,j,k,iblock)=0.0d0
346 eethet(m,l,i,j,k,iblock)=0.0d0
352 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
353 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
363 do j=-nthetyp,nthetyp
364 do k=-nthetyp,nthetyp
365 read (ithep,'(6a)') res1
366 read (ithep,*) aa0thet(i,j,k,iblock)
367 read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
369 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
370 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
371 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
372 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle)
375 & (((ffthet(llll,lll,ll,i,j,k,iblock),
376 & ffthet(lll,llll,ll,i,j,k,iblock),
377 & ggthet(llll,lll,ll,i,j,k,iblock)
378 & ,ggthet(lll,llll,ll,i,j,k,iblock),
379 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
384 C For dummy ends assign glycine-type coefficients of theta-only terms; the
385 C coefficients of theta-and-gamma-dependent terms are zero.
390 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
391 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
393 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
394 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
397 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
399 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
402 C Substitution for D aminoacids from symmetry.
405 do j=-nthetyp,nthetyp
406 do k=-nthetyp,nthetyp
407 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
409 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
413 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
414 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
415 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
416 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
422 ffthet(llll,lll,ll,i,j,k,iblock)=
423 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
424 ffthet(lll,llll,ll,i,j,k,iblock)=
425 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
426 ggthet(llll,lll,ll,i,j,k,iblock)=
427 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
428 ggthet(lll,llll,ll,i,j,k,iblock)=
429 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
439 C Control printout of the coefficients of virtual-bond-angle potentials
442 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
446 write (iout,'(//4a)')
447 & 'Type ',onelett(i),onelett(j),onelett(k)
448 write (iout,'(//a,10x,a)') " l","a[l]"
449 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
450 write (iout,'(i2,1pe15.5)')
451 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
453 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
454 & "b",l,"c",l,"d",l,"e",l
456 write (iout,'(i2,4(1pe15.5))') m,
457 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
458 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
462 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
463 & "f+",l,"f-",l,"g+",l,"g-",l
466 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
467 & ffthet(n,m,l,i,j,k,iblock),
468 & ffthet(m,n,l,i,j,k,iblock),
469 & ggthet(n,m,l,i,j,k,iblock),
470 & ggthet(m,n,l,i,j,k,iblock)
483 C Read the parameters of the probability distribution/energy expression
484 C of the side chains.
487 cc write (iout,*) "tu dochodze",i
488 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
492 dsc_inv(i)=1.0D0/dsc(i)
503 read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
504 censc(1,1,-i)=censc(1,1,i)
505 censc(2,1,-i)=censc(2,1,i)
506 censc(3,1,-i)=-censc(3,1,i)
508 read (irotam,*) bsc(j,i)
509 read (irotam,*) (censc(k,j,i),k=1,3),
510 & ((blower(k,l,j),l=1,k),k=1,3)
511 censc(1,j,-i)=censc(1,j,i)
512 censc(2,j,-i)=censc(2,j,i)
513 censc(3,j,-i)=-censc(3,j,i)
514 C BSC is amplitude of Gaussian
521 akl=akl+blower(k,m,j)*blower(l,m,j)
525 if (((k.eq.3).and.(l.ne.3))
526 & .or.((l.eq.3).and.(k.ne.3))) then
527 gaussc(k,l,j,-i)=-akl
528 gaussc(l,k,j,-i)=-akl
540 write (iout,'(/a)') 'Parameters of side-chain local geometry'
544 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
545 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
546 c write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
547 c write (iout,'(a,f10.4,4(16x,f10.4))')
548 c & 'Center ',(bsc(j,i),j=1,nlobi)
549 c write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
550 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
551 & 'log h',(bsc(j,i),j=1,nlobi)
552 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
553 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
560 c blower(k,l,j)=gaussc(ind,j,i)
565 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
566 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
573 C Read scrot parameters for potentials determined from all-atom AM1 calculations
574 C added by Urszula Kozlowska 07/11/2007
582 read(irotam,*) sc_parmin(j,i)
590 C Read torsional parameters in old format
592 read (itorp,*) ntortyp,nterm_old
593 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
594 read (itorp,*) (itortyp(i),i=1,ntyp)
599 read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
605 write (iout,'(/a/)') 'Torsional constants:'
608 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
609 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
617 C Read torsional parameters
619 read (itorp,*) ntortyp
620 read (itorp,*) (itortyp(i),i=1,ntyp)
621 write (iout,*) 'ntortyp',ntortyp
624 itortyp(i)=-itortyp(-i)
626 c write (iout,*) 'ntortyp',ntortyp
628 do j=-ntortyp+1,ntortyp-1
629 read (itorp,*) nterm(i,j,iblock),
631 nterm(-i,-j,iblock)=nterm(i,j,iblock)
632 nlor(-i,-j,iblock)=nlor(i,j,iblock)
635 do k=1,nterm(i,j,iblock)
636 read (itorp,*) kk,v1(k,i,j,iblock),
638 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
639 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
640 v0ij=v0ij+si*v1(k,i,j,iblock)
643 do k=1,nlor(i,j,iblock)
644 read (itorp,*) kk,vlor1(k,i,j),
645 & vlor2(k,i,j),vlor3(k,i,j)
646 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
649 v0(-i,-j,iblock)=v0ij
655 write (iout,'(/a/)') 'Torsional constants:'
658 write (iout,*) 'ityp',i,' jtyp',j
659 write (iout,*) 'Fourier constants'
660 do k=1,nterm(i,j,iblock)
661 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
664 write (iout,*) 'Lorenz constants'
665 do k=1,nlor(i,j,iblock)
666 write (iout,'(3(1pe15.5))')
667 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
673 C 6/23/01 Read parameters for double torsionals
677 do j=-ntortyp+1,ntortyp-1
678 do k=-ntortyp+1,ntortyp-1
679 read (itordp,'(3a1)') t1,t2,t3
680 c write (iout,*) "OK onelett",
683 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
684 & .or. t3.ne.toronelet(k)) then
685 write (iout,*) "Error in double torsional parameter file",
688 call MPI_Finalize(Ierror)
690 stop "Error in double torsional parameter file"
692 read (itordp,*) ntermd_1(i,j,k,iblock),
693 & ntermd_2(i,j,k,iblock)
694 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
695 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
696 read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,
697 & ntermd_1(i,j,k,iblock))
698 read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,
699 & ntermd_1(i,j,k,iblock))
700 read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,
701 & ntermd_1(i,j,k,iblock))
702 read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,
703 & ntermd_1(i,j,k,iblock))
704 C Martix of D parameters for one dimesional foureir series
705 do l=1,ntermd_1(i,j,k,iblock)
706 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
707 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
708 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
709 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
710 c write(iout,*) "whcodze" ,
711 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
713 read (itordp,*) ((v2c(l,m,i,j,k,iblock),
714 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
715 & v2s(m,l,i,j,k,iblock),
716 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
717 C Martix of D parameters for two dimesional fourier series
718 do l=1,ntermd_2(i,j,k,iblock)
720 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
721 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
722 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
723 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
732 write (iout,*) 'Constants for double torsionals'
735 do j=-ntortyp+1,ntortyp-1
736 do k=-ntortyp+1,ntortyp-1
737 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
738 & ' nsingle',ntermd_1(i,j,k,iblock),
739 & ' ndouble',ntermd_2(i,j,k,iblock)
741 write (iout,*) 'Single angles:'
742 do l=1,ntermd_1(i,j,k,iblock)
743 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
744 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
745 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
746 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
749 write (iout,*) 'Pairs of angles:'
750 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
751 do l=1,ntermd_2(i,j,k,iblock)
752 write (iout,'(i5,20f10.5)')
753 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
756 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
757 do l=1,ntermd_2(i,j,k,iblock)
758 write (iout,'(i5,20f10.5)')
759 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
760 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
769 C Read of Side-chain backbone correlation parameters
770 C Modified 11 May 2012 by Adasko
773 read (isccor,*) nsccortyp
774 read (isccor,*) (isccortyp(i),i=1,ntyp)
776 isccortyp(i)=-isccortyp(-i)
778 iscprol=isccortyp(20)
779 c write (iout,*) 'ntortyp',ntortyp
781 cc maxinter is maximum interaction sites
786 &nterm_sccor(i,j),nlor_sccor(i,j)
787 write (iout,*) nterm_sccor(i,j)
793 nterm_sccor(-i,j)=nterm_sccor(i,j)
794 nterm_sccor(-i,-j)=nterm_sccor(i,j)
795 nterm_sccor(i,-j)=nterm_sccor(i,j)
796 write (iout,*) nterm_sccor(i,j),nterm_sccor(-i,j),
797 & nterm_sccor(-i,-j),nterm_sccor(i,-j)
798 do k=1,nterm_sccor(i,j)
799 read (isccor,*) kk,v1sccor(k,l,i,j)
801 if (j.eq.iscprol) then
802 if (i.eq.isccortyp(10)) then
803 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
804 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
806 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
807 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
808 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
809 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
810 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
811 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
812 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
813 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
816 if (i.eq.isccortyp(10)) then
817 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
818 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
820 if (j.eq.isccortyp(10)) then
821 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
822 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
824 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
825 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
826 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
827 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
828 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
829 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
833 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
834 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
835 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
836 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
839 do k=1,nlor_sccor(i,j)
840 read (isccor,*) kk,vlor1sccor(k,i,j),
841 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
842 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
843 &(1+vlor3sccor(k,i,j)**2)
845 v0sccor(l,i,j)=v0ijsccor
846 v0sccor(l,-i,j)=v0ijsccor1
847 v0sccor(l,i,-j)=v0ijsccor2
848 v0sccor(l,-i,-j)=v0ijsccor3
854 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
857 write (iout,*) 'ityp',i,' jtyp',j
858 write (iout,*) 'Fourier constants'
859 do k=1,nterm_sccor(i,j)
860 write (iout,'(2(1pe15.5))')
861 & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
863 write (iout,*) 'Lorenz constants'
864 do k=1,nlor_sccor(i,j)
865 write (iout,'(3(1pe15.5))')
866 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
872 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
873 C interaction energy of the Gly, Ala, and Pro prototypes.
875 read (ifourier,*) nloctyp
878 read (ifourier,*) (b(ii,i),ii=1,13)
880 write (iout,*) 'Type',i
881 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
889 B1tilde(1,i) = b(3,i)
890 B1tilde(2,i) =-b(5,i)
891 B1tilde(1,-i) =-b(3,i)
892 B1tilde(2,-i) =b(5,i)
916 Ctilde(2,1,i)=-b(9,i)
918 Ctilde(1,1,-i)=b(7,i)
919 Ctilde(1,2,-i)=-b(9,i)
920 Ctilde(2,1,-i)=b(9,i)
921 Ctilde(2,2,-i)=b(7,i)
923 c Ctilde(1,1,i)=0.0d0
924 c Ctilde(1,2,i)=0.0d0
925 c Ctilde(2,1,i)=0.0d0
926 c Ctilde(2,2,i)=0.0d0
941 Dtilde(2,1,i)=-b(8,i)
943 Dtilde(1,1,-i)=b(6,i)
944 Dtilde(1,2,-i)=-b(8,i)
945 Dtilde(2,1,-i)=b(8,i)
946 Dtilde(2,2,-i)=b(6,i)
948 c Dtilde(1,1,i)=0.0d0
949 c Dtilde(1,2,i)=0.0d0
950 c Dtilde(2,1,i)=0.0d0
951 c Dtilde(2,2,i)=0.0d0
952 EE(1,1,i)= b(10,i)+b(11,i)
953 EE(2,2,i)=-b(10,i)+b(11,i)
954 EE(2,1,i)= b(12,i)-b(13,i)
955 EE(1,2,i)= b(12,i)+b(13,i)
956 EE(1,1,-i)= b(10,i)+b(11,i)
957 EE(2,2,-i)=-b(10,i)+b(11,i)
958 EE(2,1,-i)=-b(12,i)+b(13,i)
959 EE(1,2,-i)=-b(12,i)-b(13,i)
965 c ee(2,1,i)=ee(1,2,i)
970 write (iout,*) 'Type',i
972 c write (iout,'(f10.5)') B1(:,i)
973 write(iout,*) B1(1,i),B1(2,i)
975 c write (iout,'(f10.5)') B2(:,i)
976 write(iout,*) B2(1,i),B2(2,i)
979 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
983 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
987 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
992 C Read electrostatic-interaction parameters
995 write (iout,'(/a)') 'Electrostatic interaction constants:'
996 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
997 & 'IT','JT','APP','BPP','AEL6','AEL3'
999 read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
1000 read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
1001 read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
1002 read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
1007 app (i,j)=epp(i,j)*rri*rri
1008 bpp (i,j)=-2.0D0*epp(i,j)*rri
1009 ael6(i,j)=elpp6(i,j)*4.2D0**6
1010 ael3(i,j)=elpp3(i,j)*4.2D0**3
1012 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1013 & ael6(i,j),ael3(i,j)
1018 C Read side-chain interaction parameters.
1020 read (isidep,*) ipot,expon
1021 if (ipot.lt.1 .or. ipot.gt.5) then
1022 write (iout,'(2a)') 'Error while reading SC interaction',
1023 & 'potential file - unknown potential type.'
1027 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1028 & ', exponents are ',expon,2*expon
1029 goto (10,20,30,30,40) ipot
1030 C----------------------- LJ potential ---------------------------------
1031 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1033 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1034 write (iout,'(a/)') 'The epsilon array:'
1035 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1036 write (iout,'(/a)') 'One-body parameters:'
1037 write (iout,'(a,4x,a)') 'residue','sigma'
1038 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1041 C----------------------- LJK potential --------------------------------
1042 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1043 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1045 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1046 write (iout,'(a/)') 'The epsilon array:'
1047 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1048 write (iout,'(/a)') 'One-body parameters:'
1049 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1050 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1054 C---------------------- GB or BP potential -----------------------------
1056 read (isidep,*)(eps(i,j),j=i,ntyp)
1058 read (isidep,*)(sigma0(i),i=1,ntyp)
1059 read (isidep,*)(sigii(i),i=1,ntyp)
1060 read (isidep,*)(chip(i),i=1,ntyp)
1061 read (isidep,*)(alp(i),i=1,ntyp)
1062 C For the GB potential convert sigma'**2 into chi'
1065 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1069 write (iout,'(/a/)') 'Parameters of the BP potential:'
1070 write (iout,'(a/)') 'The epsilon array:'
1071 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1072 write (iout,'(/a)') 'One-body parameters:'
1073 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1075 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1076 & chip(i),alp(i),i=1,ntyp)
1079 C--------------------- GBV potential -----------------------------------
1080 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1081 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1082 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1084 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1085 write (iout,'(a/)') 'The epsilon array:'
1086 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1087 write (iout,'(/a)') 'One-body parameters:'
1088 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1089 & 's||/s_|_^2',' chip ',' alph '
1090 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1091 & sigii(i),chip(i),alp(i),i=1,ntyp)
1095 C-----------------------------------------------------------------------
1096 C Calculate the "working" parameters of SC interactions.
1104 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1105 sigma(j,i)=sigma(i,j)
1106 rs0(i,j)=dwa16*sigma(i,j)
1110 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1111 & 'Working parameters of the SC interactions:',
1112 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1117 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1126 sigeps=dsign(1.0D0,epsij)
1128 aa(i,j)=epsij*rrij*rrij
1129 bb(i,j)=-sigeps*epsij*rrij
1133 sigt1sq=sigma0(i)**2
1134 sigt2sq=sigma0(j)**2
1137 ratsig1=sigt2sq/sigt1sq
1138 ratsig2=1.0D0/ratsig1
1139 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1140 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1141 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1145 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1146 sigmaii(i,j)=rsum_max
1147 sigmaii(j,i)=rsum_max
1149 c sigmaii(i,j)=r0(i,j)
1150 c sigmaii(j,i)=r0(i,j)
1152 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1153 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1154 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1155 augm(i,j)=epsij*r_augm**(2*expon)
1156 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1163 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1164 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1165 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1170 C Define the SC-p interaction constants
1174 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1176 c aad(i,1)=0.3D0*4.0D0**12
1177 C Following line for constants currently implemented
1178 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1179 aad(i,1)=1.5D0*4.0D0**12
1180 c aad(i,1)=0.17D0*5.6D0**12
1182 C "Soft" SC-p repulsion
1184 C Following line for constants currently implemented
1185 c aad(i,1)=0.3D0*4.0D0**6
1186 C "Hard" SC-p repulsion
1187 bad(i,1)=3.0D0*4.0D0**6
1188 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1197 C 8/9/01 Read the SC-p interaction constants from file
1200 read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1203 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1204 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1205 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1206 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1210 write (iout,*) "Parameters of SC-p interactions:"
1212 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1213 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1218 C Define the constants of the disulfide bridge
1222 c Old arbitrary potential - commented out.
1227 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1228 c energy surface of diethyl disulfide.
1229 c A. Liwo and U. Kozlowska, 11/24/03
1240 write (iout,'(/a)') "Disulfide bridge parameters:"
1241 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1242 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1243 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1244 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,