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"/
27 dimension blower(3,3,maxlob)
28 character*800 controlcard
29 character*256 bondname_t,thetname_t,rotname_t,torname_t,
30 & tordname_t,fouriername_t,elename_t,sidename_t,scpname_t,
36 double precision ip,mp
40 C Set LPRINT=.TRUE. for debugging
41 dwa16=2.0d0**(1.0d0/6.0d0)
44 C Assign virtual-bond length
48 call card_concat(controlcard,.true.)
51 key = wname(i)(:ilen(wname(i)))
52 call reada(controlcard,key(:ilen(key)),ww(i),1.0d0)
54 call reada(controlcard,"D0CM",d0cm,3.78d0)
55 call reada(controlcard,"AKCM",akcm,15.1d0)
56 call reada(controlcard,"AKTH",akth,11.0d0)
57 call reada(controlcard,"AKCT",akct,12.0d0)
58 call reada(controlcard,"V1SS",v1ss,-1.08d0)
59 call reada(controlcard,"V2SS",v2ss,7.61d0)
60 call reada(controlcard,"V3SS",v3ss,13.7d0)
61 call reada(controlcard,"EBR",ebr,-5.50D0)
62 c dyn_ss=(index(controlcard,'DYN_SS').gt.0)
63 write (iout,*) "iparm",iparm," myparm",myparm
65 c dyn_ss_mask(i)=.false.
69 dyn_ssbond_ij(i,j)=1.0d300
72 call reada(controlcard,"HT",Ht,0.0D0)
74 c if(me.eq.king.or..not.out1file) then
75 c print *,'indpdb=',indpdb,' pdbref=',pdbref
77 c If reading not own parameters, skip assignment
78 cc write(iout,*) "KURWA", ww(15)
79 if (iparm.eq.myparm .or. .not.separate_parset) then
82 c Setup weights for UNRES
104 cc write(iout,*) "KURWA", wstrain,akcm,akth,wsc,dyn_ss
106 call card_concat(controlcard,.false.)
108 c Return if not own parameters
110 if (iparm.ne.myparm .and. separate_parset) return
112 call reads(controlcard,"BONDPAR",bondname_t,bondname)
113 open (ibond,file=bondname_t,status='old')
115 call reads(controlcard,"THETPAR",thetname_t,thetname)
116 open (ithep,file=thetname_t,status='old')
118 call reads(controlcard,"ROTPAR",rotname_t,rotname)
119 open (irotam,file=rotname_t,status='old')
121 call reads(controlcard,"TORPAR",torname_t,torname)
122 open (itorp,file=torname_t,status='old')
124 call reads(controlcard,"TORDPAR",tordname_t,tordname)
125 open (itordp,file=tordname_t,status='old')
127 call reads(controlcard,"SCCORPAR",sccorname_t,sccorname)
128 open (isccor,file=sccorname_t,status='old')
130 call reads(controlcard,"FOURIER",fouriername_t,fouriername)
131 open (ifourier,file=fouriername_t,status='old')
133 call reads(controlcard,"ELEPAR",elename_t,elename)
134 open (ielep,file=elename_t,status='old')
136 call reads(controlcard,"SIDEPAR",sidename_t,sidename)
137 open (isidep,file=sidename_t,status='old')
139 call reads(controlcard,"SCPPAR",scpname_t,scpname)
140 open (iscpp,file=scpname_t,status='old')
142 write (iout,*) "Parameter set:",iparm
143 write (iout,*) "Energy-term weights:"
145 write (iout,'(a16,f10.5)') wname(i),ww(i)
147 write (iout,*) "Sidechain potential file : ",
148 & sidename_t(:ilen(sidename_t))
150 write (iout,*) "SCp potential file : ",
151 & scpname_t(:ilen(scpname_t))
153 write (iout,*) "Electrostatic potential file : ",
154 & elename_t(:ilen(elename_t))
155 write (iout,*) "Cumulant coefficient file : ",
156 & fouriername_t(:ilen(fouriername_t))
157 write (iout,*) "Torsional parameter file : ",
158 & torname_t(:ilen(torname_t))
159 write (iout,*) "Double torsional parameter file : ",
160 & tordname_t(:ilen(tordname_t))
161 write (iout,*) "Backbone-rotamer parameter file : ",
162 & sccorname(:ilen(sccorname))
163 write (iout,*) "Bond & inertia constant file : ",
164 & bondname_t(:ilen(bondname_t))
165 write (iout,*) "Bending parameter file : ",
166 & thetname_t(:ilen(thetname_t))
167 write (iout,*) "Rotamer parameter file : ",
168 & rotname_t(:ilen(rotname_t))
171 c Read the virtual-bond parameters, masses, and moments of inertia
172 c and Stokes' radii of the peptide group and side chains
175 read (ibond,*) vbldp0,akp
178 read (ibond,*) vbldsc0(1,i),aksc(1,i)
179 dsc(i) = vbldsc0(1,i)
183 dsc_inv(i)=1.0D0/dsc(i)
187 read (ibond,*) ijunk,vbldp0,akp,rjunk
189 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
191 dsc(i) = vbldsc0(1,i)
195 dsc_inv(i)=1.0D0/dsc(i)
200 write(iout,'(/a/)')"Force constants virtual bonds:"
201 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
203 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
205 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
206 & vbldsc0(1,i),aksc(1,i),abond0(1,i)
208 write (iout,'(13x,3f10.5)')
209 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
215 C Read the parameters of the probability distribution/energy expression
216 C of the virtual-bond valence angles theta
219 read (ithep,*) a0thet(i),(athet(j,i,1,1),j=1,2),
220 & (bthet(j,i,1,1),j=1,2)
221 read (ithep,*) (polthet(j,i),j=0,3)
222 read (ithep,*) (gthet(j,i),j=1,3)
223 read (ithep,*) theta0(i),sig0(i),sigc0(i)
227 athet(1,i,1,-1)=athet(1,i,1,1)
228 athet(2,i,1,-1)=athet(2,i,1,1)
229 bthet(1,i,1,-1)=-bthet(1,i,1,1)
230 bthet(2,i,1,-1)=-bthet(2,i,1,1)
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)
238 athet(1,i,-1,-1)=athet(1,-i,1,1)
239 athet(2,i,-1,-1)=-athet(2,-i,1,1)
240 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
241 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)
254 polthet(j,i)=polthet(j,-i)
257 gthet(j,i)=gthet(j,-i)
263 c & 'Parameters of the virtual-bond valence angles:'
264 c write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
265 c & ' ATHETA0 ',' A1 ',' A2 ',
268 c write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
269 c & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
271 c write (iout,'(/a/9x,5a/79(1h-))')
272 c & 'Parameters of the expression for sigma(theta_c):',
273 c & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
274 c & ' ALPH3 ',' SIGMA0C '
276 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
277 c & (polthet(j,i),j=0,3),sigc0(i)
279 c write (iout,'(/a/9x,5a/79(1h-))')
280 c & 'Parameters of the second gaussian:',
281 c & ' THETA0 ',' SIGMA0 ',' G1 ',
284 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
285 c & sig0(i),(gthet(j,i),j=1,3)
288 & 'Parameters of the virtual-bond valence angles:'
289 write (iout,'(/a/9x,5a/79(1h-))')
290 & 'Coefficients of expansion',
291 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
292 & ' b1*10^1 ',' b2*10^1 '
294 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
295 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
296 & (10*bthet(j,i,1,1),j=1,2)
298 write (iout,'(/a/9x,5a/79(1h-))')
299 & 'Parameters of the expression for sigma(theta_c):',
300 & ' alpha0 ',' alph1 ',' alph2 ',
301 & ' alhp3 ',' sigma0c '
303 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
304 & (polthet(j,i),j=0,3),sigc0(i)
306 write (iout,'(/a/9x,5a/79(1h-))')
307 & 'Parameters of the second gaussian:',
308 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
311 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
312 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
317 C Read the parameters of Utheta determined from ab initio surfaces
318 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
320 read (ithep,*) nthetyp,ntheterm,ntheterm2,
321 & ntheterm3,nsingle,ndouble
322 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
323 read (ithep,*) (ithetyp(i),i=1,ntyp1)
325 ithetyp(i)=-ithetyp(-i)
327 c write (iout,*) "tu dochodze"
329 do i=-maxthetyp,maxthetyp
330 do j=-maxthetyp,maxthetyp
331 do k=-maxthetyp,maxthetyp
332 aa0thet(i,j,k,iblock)=0.0d0
334 aathet(l,i,j,k,iblock)=0.0d0
338 bbthet(m,l,i,j,k,iblock)=0.0d0
339 ccthet(m,l,i,j,k,iblock)=0.0d0
340 ddthet(m,l,i,j,k,iblock)=0.0d0
341 eethet(m,l,i,j,k,iblock)=0.0d0
347 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
348 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
359 do j=-nthetyp,nthetyp
360 do k=-nthetyp,nthetyp
361 read (ithep,'(6a)') res1
362 read (ithep,*) aa0thet(i,j,k,iblock)
363 read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
365 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
366 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
367 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
368 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle)
371 & (((ffthet(llll,lll,ll,i,j,k,iblock),
372 & ffthet(lll,llll,ll,i,j,k,iblock),
373 & ggthet(llll,lll,ll,i,j,k,iblock)
374 & ,ggthet(lll,llll,ll,i,j,k,iblock),
375 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
380 C For dummy ends assign glycine-type coefficients of theta-only terms; the
381 C coefficients of theta-and-gamma-dependent terms are zero.
386 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
387 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
389 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
390 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
393 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
395 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
398 C Substitution for D aminoacids from symmetry.
401 do j=-nthetyp,nthetyp
402 do k=-nthetyp,nthetyp
403 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
405 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
409 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
410 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
411 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
412 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
418 ffthet(llll,lll,ll,i,j,k,iblock)=
419 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
420 ffthet(lll,llll,ll,i,j,k,iblock)=
421 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
422 ggthet(llll,lll,ll,i,j,k,iblock)=
423 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
424 ggthet(lll,llll,ll,i,j,k,iblock)=
425 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
434 C Control printout of the coefficients of virtual-bond-angle potentials
437 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
441 write (iout,'(//4a)')
442 & 'Type ',onelett(i),onelett(j),onelett(k)
443 write (iout,'(//a,10x,a)') " l","a[l]"
444 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
445 write (iout,'(i2,1pe15.5)')
446 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
448 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
449 & "b",l,"c",l,"d",l,"e",l
451 write (iout,'(i2,4(1pe15.5))') m,
452 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
453 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
457 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
458 & "f+",l,"f-",l,"g+",l,"g-",l
461 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
462 & ffthet(n,m,l,i,j,k,iblock),
463 & ffthet(m,n,l,i,j,k,iblock),
464 & ggthet(n,m,l,i,j,k,iblock),
465 & ggthet(m,n,l,i,j,k,iblock)
478 C Read the parameters of the probability distribution/energy expression
479 C of the side chains.
482 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
486 dsc_inv(i)=1.0D0/dsc(i)
497 read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
498 censc(1,1,-i)=censc(1,1,i)
499 censc(2,1,-i)=censc(2,1,i)
500 censc(3,1,-i)=-censc(3,1,i)
502 read (irotam,*) bsc(j,i)
503 read (irotam,*) (censc(k,j,i),k=1,3),
504 & ((blower(k,l,j),l=1,k),k=1,3)
505 censc(1,j,-i)=censc(1,j,i)
506 censc(2,j,-i)=censc(2,j,i)
507 censc(3,j,-i)=-censc(3,j,i)
514 akl=akl+blower(k,m,j)*blower(l,m,j)
518 if (((k.eq.3).and.(l.ne.3))
519 & .or.((l.eq.3).and.(k.ne.3))) then
520 gaussc(k,l,j,-i)=-akl
521 gaussc(l,k,j,-i)=-akl
533 write (iout,'(/a)') 'Parameters of side-chain local geometry'
537 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
538 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
539 c write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
540 c write (iout,'(a,f10.4,4(16x,f10.4))')
541 c & 'Center ',(bsc(j,i),j=1,nlobi)
542 c write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
543 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
544 & 'log h',(bsc(j,i),j=1,nlobi)
545 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
546 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
553 c blower(k,l,j)=gaussc(ind,j,i)
558 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
559 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
566 C Read scrot parameters for potentials determined from all-atom AM1 calculations
567 C added by Urszula Kozlowska 07/11/2007
575 read(irotam,*) sc_parmin(j,i)
583 C Read torsional parameters in old format
585 read (itorp,*) ntortyp,nterm_old
586 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
587 read (itorp,*) (itortyp(i),i=1,ntyp)
592 read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
598 write (iout,'(/a/)') 'Torsional constants:'
601 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
602 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
610 C Read torsional parameters
612 read (itorp,*) ntortyp
613 read (itorp,*) (itortyp(i),i=1,ntyp)
616 itortyp(i)=-itortyp(-i)
618 c write (iout,*) 'ntortyp',ntortyp
620 do j=-ntortyp+1,ntortyp-1
621 read (itorp,*) nterm(i,j,iblock),
623 nterm(-i,-j,iblock)=nterm(i,j,iblock)
624 nlor(-i,-j,iblock)=nlor(i,j,iblock)
627 do k=1,nterm(i,j,iblock)
628 read (itorp,*) kk,v1(k,i,j,iblock),v2(k,i,j,iblock)
629 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
630 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
631 v0ij=v0ij+si*v1(k,i,j,iblock)
634 do k=1,nlor(i,j,iblock)
635 read (itorp,*) kk,vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
636 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
639 v0(-i,-j,iblock)=v0ij
645 write (iout,'(/a/)') 'Torsional constants:'
648 write (iout,*) 'ityp',i,' jtyp',j
649 write (iout,*) 'Fourier constants'
650 do k=1,nterm(i,j,iblock)
651 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
654 write (iout,*) 'Lorenz constants'
655 do k=1,nlor(i,j,iblock)
656 write (iout,'(3(1pe15.5))')
657 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
663 C 6/23/01 Read parameters for double torsionals
667 do j=-ntortyp+1,ntortyp-1
668 do k=-ntortyp+1,ntortyp-1
669 read (itordp,'(3a1)') t1,t2,t3
670 if (t1.ne.onelett(i) .or. t2.ne.onelett(j)
671 & .or. t3.ne.onelett(k)) then
672 write (iout,*) "Error in double torsional parameter file",
674 stop "Error in double torsional parameter file"
676 read (itordp,*) ntermd_1(i,j,k,iblock),
677 & ntermd_2(i,j,k,iblock)
678 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
679 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
680 read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,
681 & ntermd_1(i,j,k,iblock))
682 read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,
683 & ntermd_1(i,j,k,iblock))
684 read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,
685 & ntermd_1(i,j,k,iblock))
686 read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,
687 & ntermd_1(i,j,k,iblock))
688 C Martix of D parameters for one dimesional foureir series
689 do l=1,ntermd_1(i,j,k,iblock)
690 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
691 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
692 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
693 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
694 c write(iout,*) "whcodze" ,
695 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
697 read (itordp,*) ((v2c(l,m,i,j,k,iblock),
698 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
699 & v2s(m,l,i,j,k,iblock),
700 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
701 C Martix of D parameters for two dimesional fourier series
702 do l=1,ntermd_2(i,j,k,iblock)
704 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
705 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
706 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
707 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
716 write (iout,*) 'Constants for double torsionals'
719 do j=-ntortyp+1,ntortyp-1
720 do k=-ntortyp+1,ntortyp-1
721 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
722 & ' nsingle',ntermd_1(i,j,k,iblock),
723 & ' ndouble',ntermd_2(i,j,k,iblock)
725 write (iout,*) 'Single angles:'
726 do l=1,ntermd_1(i,j,k,iblock)
727 write (iout,'(i5,2f10.5,5x,2f10.5)') l,
728 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
729 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock)
732 write (iout,*) 'Pairs of angles:'
733 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
734 do l=1,ntermd_2(i,j,k,iblock)
735 write (iout,'(i5,20f10.5)')
736 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
739 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
740 do l=1,ntermd_2(i,j,k,iblock)
741 write (iout,'(i5,20f10.5)')
742 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
751 C Read of Side-chain backbone correlation parameters
752 C Modified 11 May 2012 by Adasko
755 read (isccor,*) nsccortyp
756 read (isccor,*) (isccortyp(i),i=1,ntyp)
758 isccortyp(i)=-isccortyp(-i)
760 iscprol=isccortyp(20)
761 c write (iout,*) 'ntortyp',ntortyp
763 cc maxinter is maximum interaction sites
767 read (isccor,*) nterm_sccor(i,j),nlor_sccor(i,j)
773 nterm_sccor(-i,j)=nterm_sccor(i,j)
774 nterm_sccor(-i,-j)=nterm_sccor(i,j)
775 nterm_sccor(i,-j)=nterm_sccor(i,j)
776 do k=1,nterm_sccor(i,j)
777 read (isccor,*) kk,v1sccor(k,l,i,j)
779 if (j.eq.iscprol) then
780 if (i.eq.isccortyp(10)) then
781 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
782 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
784 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
785 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
786 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
787 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
788 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
789 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
790 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
791 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
794 if (i.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 if (j.eq.isccortyp(10)) then
799 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
800 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)
804 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
805 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
806 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
807 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
811 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
812 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
813 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
814 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
817 do k=1,nlor_sccor(i,j)
818 read (isccor,*) kk,vlor1sccor(k,i,j),
819 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
820 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
821 &(1+vlor3sccor(k,i,j)**2)
823 v0sccor(l,i,j)=v0ijsccor
824 v0sccor(l,-i,j)=v0ijsccor1
825 v0sccor(l,i,-j)=v0ijsccor2
826 v0sccor(l,-i,-j)=v0ijsccor3
833 write (iout,'(/a/)') 'Torsional constants:'
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))')v1sccor(k,l,i,j),v2sccor(k,l,i,j)
841 write (iout,*) 'Lorenz constants'
842 do k=1,nlor_sccor(i,j)
843 write (iout,'(3(1pe15.5))')
844 & 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),ii=1,13)
859 write (iout,*) 'Type',i
860 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
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
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)+b(11)
932 EE(2,2,i)=-b(10)+b(11)
933 EE(2,1,i)= b(12)-b(13)
934 EE(1,2,i)= b(12)+b(13)
935 EE(1,1,-i)= b(10)+b(11)
936 EE(2,2,-i)=-b(10)+b(11)
937 EE(2,1,-i)=-b(12)+b(13)
938 EE(1,2,-i)=-b(12)-b(13)
944 c ee(2,1,i)=ee(1,2,i)
948 write (iout,*) 'Type',i
950 c write (iout,'(f10.5)') B1(:,i)
951 write(iout,*) B1(1,i),B1(2,i)
953 c write (iout,'(f10.5)') B2(:,i)
954 write(iout,*) B2(1,i),B2(2,i)
957 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
961 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
965 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
970 C Read electrostatic-interaction parameters
973 write (iout,'(/a)') 'Electrostatic interaction constants:'
974 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
975 & 'IT','JT','APP','BPP','AEL6','AEL3'
977 read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
978 read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
979 read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
980 read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
985 app (i,j)=epp(i,j)*rri*rri
986 bpp (i,j)=-2.0D0*epp(i,j)*rri
987 ael6(i,j)=elpp6(i,j)*4.2D0**6
988 ael3(i,j)=elpp3(i,j)*4.2D0**3
989 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
990 & ael6(i,j),ael3(i,j)
994 C Read side-chain interaction parameters.
996 read (isidep,*) ipot,expon
997 if (ipot.lt.1 .or. ipot.gt.5) then
998 write (iout,'(2a)') 'Error while reading SC interaction',
999 & 'potential file - unknown potential type.'
1003 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1004 & ', exponents are ',expon,2*expon
1005 goto (10,20,30,30,40) ipot
1006 C----------------------- LJ potential ---------------------------------
1007 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1009 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1010 write (iout,'(a/)') 'The epsilon array:'
1011 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1012 write (iout,'(/a)') 'One-body parameters:'
1013 write (iout,'(a,4x,a)') 'residue','sigma'
1014 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1017 C----------------------- LJK potential --------------------------------
1018 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1019 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1021 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1022 write (iout,'(a/)') 'The epsilon array:'
1023 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1024 write (iout,'(/a)') 'One-body parameters:'
1025 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1026 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1030 C---------------------- GB or BP potential -----------------------------
1031 30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1032 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(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 and SS bond potentials
1145 ss_depth=ebr/wsc-0.25*eps(1,1)
1146 Ht=Ht/wsc-0.25*eps(1,1)
1147 akcm=akcm*wstrain/wsc
1148 akth=akth*wstrain/wsc
1149 akct=akct*wstrain/wsc
1150 v1ss=v1ss*wstrain/wsc
1151 v2ss=v2ss*wstrain/wsc
1152 v3ss=v3ss*wstrain/wsc
1154 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
1156 write (iout,*) "Parameters of the SS-bond potential:"
1157 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
1159 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
1160 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
1161 write (iout,*)" HT",Ht
1165 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1167 c aad(i,1)=0.3D0*4.0D0**12
1168 C Following line for constants currently implemented
1169 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1170 aad(i,1)=1.5D0*4.0D0**12
1171 c aad(i,1)=0.17D0*5.6D0**12
1173 C "Soft" SC-p repulsion
1175 C Following line for constants currently implemented
1176 c aad(i,1)=0.3D0*4.0D0**6
1177 C "Hard" SC-p repulsion
1178 bad(i,1)=3.0D0*4.0D0**6
1179 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1188 C 8/9/01 Read the SC-p interaction constants from file
1191 read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1194 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1195 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1196 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1197 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1201 write (iout,*) "Parameters of SC-p interactions:"
1203 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1204 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1209 C Define the constants of the disulfide bridge
1213 c Old arbitrary potential - commented out.
1218 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1219 c energy surface of diethyl disulfide.
1220 c A. Liwo and U. Kozlowska, 11/24/03
1231 c write (iout,'(/a)') "Disulfide bridge parameters:"
1232 c write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1233 c write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1234 c write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1235 c write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,