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 c 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 c 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'
421 write (iout,'(//4a)')
422 & 'Type ',onelett(i),onelett(j),onelett(k)
423 write (iout,'(//a,10x,a)') " l","a[l]"
424 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
425 write (iout,'(i2,1pe15.5)')
426 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
428 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
429 & "b",l,"c",l,"d",l,"e",l
431 write (iout,'(i2,4(1pe15.5))') m,
432 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
433 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
437 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
438 & "f+",l,"f-",l,"g+",l,"g-",l
441 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
442 & ffthet(n,m,l,i,j,k,iblock),
443 & ffthet(m,n,l,i,j,k,iblock),
444 & ggthet(n,m,l,i,j,k,iblock),
445 & ggthet(m,n,l,i,j,k,iblock)
459 C Read the parameters of the probability distribution/energy expression
460 C of the side chains.
463 cc write (iout,*) "tu dochodze",i
464 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
468 dsc_inv(i)=1.0D0/dsc(i)
479 read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
480 censc(1,1,-i)=censc(1,1,i)
481 censc(2,1,-i)=censc(2,1,i)
482 censc(3,1,-i)=-censc(3,1,i)
484 read (irotam,*) bsc(j,i)
485 read (irotam,*) (censc(k,j,i),k=1,3),
486 & ((blower(k,l,j),l=1,k),k=1,3)
487 censc(1,j,-i)=censc(1,j,i)
488 censc(2,j,-i)=censc(2,j,i)
489 censc(3,j,-i)=-censc(3,j,i)
490 C BSC is amplitude of Gaussian
497 akl=akl+blower(k,m,j)*blower(l,m,j)
501 if (((k.eq.3).and.(l.ne.3))
502 & .or.((l.eq.3).and.(k.ne.3))) then
503 gaussc(k,l,j,-i)=-akl
504 gaussc(l,k,j,-i)=-akl
516 write (iout,'(/a)') 'Parameters of side-chain local geometry'
520 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
521 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
522 c write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
523 c write (iout,'(a,f10.4,4(16x,f10.4))')
524 c & 'Center ',(bsc(j,i),j=1,nlobi)
525 c write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
526 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
527 & 'log h',(bsc(j,i),j=1,nlobi)
528 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
529 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
536 c blower(k,l,j)=gaussc(ind,j,i)
541 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
542 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
549 C Read scrot parameters for potentials determined from all-atom AM1 calculations
550 C added by Urszula Kozlowska 07/11/2007
558 read(irotam,*) sc_parmin(j,i)
566 C Read torsional parameters in old format
568 read (itorp,*) ntortyp,nterm_old
569 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
570 read (itorp,*) (itortyp(i),i=1,ntyp)
575 read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
581 write (iout,'(/a/)') 'Torsional constants:'
584 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
585 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
593 C Read torsional parameters
595 read (itorp,*) ntortyp
596 read (itorp,*) (itortyp(i),i=1,ntyp)
597 write (iout,*) 'ntortyp',ntortyp
600 itortyp(i)=-itortyp(-i)
602 c write (iout,*) 'ntortyp',ntortyp
604 do j=-ntortyp+1,ntortyp-1
605 read (itorp,*) nterm(i,j,iblock),
607 nterm(-i,-j,iblock)=nterm(i,j,iblock)
608 nlor(-i,-j,iblock)=nlor(i,j,iblock)
611 do k=1,nterm(i,j,iblock)
612 read (itorp,*) kk,v1(k,i,j,iblock),
614 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
615 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
616 v0ij=v0ij+si*v1(k,i,j,iblock)
619 do k=1,nlor(i,j,iblock)
620 read (itorp,*) kk,vlor1(k,i,j),
621 & vlor2(k,i,j),vlor3(k,i,j)
622 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
625 v0(-i,-j,iblock)=v0ij
631 write (iout,'(/a/)') 'Torsional constants:'
634 write (iout,*) 'ityp',i,' jtyp',j
635 write (iout,*) 'Fourier constants'
637 do k=1,nterm(i,j,iblock)
638 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
641 write (iout,*) 'Lorenz constants'
642 do k=1,nlor(i,j,iblock)
643 write (iout,'(3(1pe15.5))')
644 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
651 C 6/23/01 Read parameters for double torsionals
655 do j=-ntortyp+1,ntortyp-1
656 do k=-ntortyp+1,ntortyp-1
657 read (itordp,'(3a1)') t1,t2,t3
658 c write (iout,*) "OK onelett",
661 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
662 & .or. t3.ne.toronelet(k)) then
663 write (iout,*) "Error in double torsional parameter file",
666 call MPI_Finalize(Ierror)
668 stop "Error in double torsional parameter file"
670 read (itordp,*) ntermd_1(i,j,k,iblock),
671 & ntermd_2(i,j,k,iblock)
672 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
673 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
674 read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,
675 & ntermd_1(i,j,k,iblock))
676 read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,
677 & ntermd_1(i,j,k,iblock))
678 read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,
679 & ntermd_1(i,j,k,iblock))
680 read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,
681 & ntermd_1(i,j,k,iblock))
682 C Martix of D parameters for one dimesional foureir series
683 do l=1,ntermd_1(i,j,k,iblock)
684 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
685 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
686 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
687 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
688 c write(iout,*) "whcodze" ,
689 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
691 read (itordp,*) ((v2c(l,m,i,j,k,iblock),
692 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
693 & v2s(m,l,i,j,k,iblock),
694 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
695 C Martix of D parameters for two dimesional fourier series
696 do l=1,ntermd_2(i,j,k,iblock)
698 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
699 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
700 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
701 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
710 write (iout,*) 'Constants for double torsionals'
713 do j=-ntortyp+1,ntortyp-1
714 do k=-ntortyp+1,ntortyp-1
715 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
716 & ' nsingle',ntermd_1(i,j,k,iblock),
717 & ' ndouble',ntermd_2(i,j,k,iblock)
719 write (iout,*) 'Single angles:'
720 do l=1,ntermd_1(i,j,k,iblock)
721 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
722 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
723 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
724 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
727 write (iout,*) 'Pairs of angles:'
728 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
729 do l=1,ntermd_2(i,j,k,iblock)
730 write (iout,'(i5,20f10.5)')
731 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
734 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
735 do l=1,ntermd_2(i,j,k,iblock)
736 write (iout,'(i5,20f10.5)')
737 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
738 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
747 C Read of Side-chain backbone correlation parameters
748 C Modified 11 May 2012 by Adasko
751 read (isccor,*) nsccortyp
752 read (isccor,*) (isccortyp(i),i=1,ntyp)
754 isccortyp(i)=-isccortyp(-i)
756 iscprol=isccortyp(20)
757 c write (iout,*) 'ntortyp',ntortyp
759 cc maxinter is maximum interaction sites
764 &nterm_sccor(i,j),nlor_sccor(i,j)
765 write (iout,*) nterm_sccor(i,j)
771 nterm_sccor(-i,j)=nterm_sccor(i,j)
772 nterm_sccor(-i,-j)=nterm_sccor(i,j)
773 nterm_sccor(i,-j)=nterm_sccor(i,j)
774 write (iout,*) 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
832 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
835 write (iout,*) 'ityp',i,' jtyp',j
836 write (iout,*) 'Fourier constants'
838 do k=1,nterm_sccor(i,j)
839 write (iout,'(2(1pe15.5))')
840 & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
843 write (iout,*) 'Lorenz constants'
844 do k=1,nlor_sccor(i,j)
845 write (iout,'(3(1pe15.5))')
846 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
852 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
853 C interaction energy of the Gly, Ala, and Pro prototypes.
855 write (iout,*) "Reading Fourier"
857 read (ifourier,*) nloctyp
860 read (ifourier,*) (b(ii,i),ii=1,13)
862 write (iout,*) 'Type',i
863 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
865 write (iout,*) "Reading Fourier 1"
868 read (ifourier,*) (bnew1(ii,1,i),ii=1,3)
869 read (ifourier,*) (bnew2(ii,1,i),ii=1,3)
870 read (ifourier,*) (bnew1(ii,2,i),ii=1,1)
871 read (ifourier,*) (bnew2(ii,2,i),ii=1,1)
872 read (ifourier,*) (eenew(ii,i),ii=1,1)
883 Ctilde(2,1,i)=-b(9,i)
885 Ctilde(1,1,-i)=b(7,i)
886 Ctilde(1,2,-i)=-b(9,i)
887 Ctilde(2,1,-i)=b(9,i)
888 Ctilde(2,2,-i)=b(7,i)
900 Dtilde(2,1,i)=-b(8,i)
902 Dtilde(1,1,-i)=b(6,i)
903 Dtilde(1,2,-i)=-b(8,i)
904 Dtilde(2,1,-i)=b(8,i)
905 Dtilde(2,2,-i)=b(6,i)
907 EEold(1,1,i)= b(10,i)+b(11,i)
908 eeold(1,1,i)= eenew(1,i)
909 EEold(2,2,i)=-b(10,i)+b(11,i)
910 EEold(2,1,i)= b(12,i)-b(13,i)
911 EEold(1,2,i)= b(12,i)+b(13,i)
912 EEold(1,1,-i)= b(10,i)+b(11,i)
913 EEold(2,2,-i)=-b(10,i)+b(11,i)
914 EEold(2,1,-i)=-b(12,i)+b(13,i)
915 EEold(1,2,-i)=-b(12,i)-b(13,i)
916 c write(iout,*) "TU DOCHODZE"
921 write (iout,*) 'Type',i
922 write (iout,*) 'BNEW1(1)'
923 write (iout,*) (BNEW1(ii,1,i),ii=1,3)
924 write (iout,*) 'BNEW1(2)'
925 write (iout,*) BNEW1(1,2,i)
926 write (iout,*) 'BNEW2(1)'
927 write (iout,*) (BNEW2(ii,1,i),ii=1,3)
928 write (iout,*) 'BNEW2(2)'
929 write (iout,*) BNEW2(1,2,i)
932 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
936 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
940 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
951 B1tilde(1,i) = b(3,i)
952 B1tilde(2,i) =-b(5,i)
953 B1tilde(1,-i) =-b(3,i)
954 B1tilde(2,-i) =b(5,i)
978 Ctilde(2,1,i)=-b(9,i)
980 Ctilde(1,1,-i)=b(7,i)
981 Ctilde(1,2,-i)=-b(9,i)
982 Ctilde(2,1,-i)=b(9,i)
983 Ctilde(2,2,-i)=b(7,i)
985 c Ctilde(1,1,i)=0.0d0
986 c Ctilde(1,2,i)=0.0d0
987 c Ctilde(2,1,i)=0.0d0
988 c Ctilde(2,2,i)=0.0d0
1001 Dtilde(1,1,i)=b(6,i)
1002 Dtilde(1,2,i)=b(8,i)
1003 Dtilde(2,1,i)=-b(8,i)
1004 Dtilde(2,2,i)=b(6,i)
1005 Dtilde(1,1,-i)=b(6,i)
1006 Dtilde(1,2,-i)=-b(8,i)
1007 Dtilde(2,1,-i)=b(8,i)
1008 Dtilde(2,2,-i)=b(6,i)
1010 c Dtilde(1,1,i)=0.0d0
1011 c Dtilde(1,2,i)=0.0d0
1012 c Dtilde(2,1,i)=0.0d0
1013 c Dtilde(2,2,i)=0.0d0
1014 EE(1,1,i)= b(10,i)+b(11,i)
1015 EE(2,2,i)=-b(10,i)+b(11,i)
1016 EE(2,1,i)= b(12,i)-b(13,i)
1017 EE(1,2,i)= b(12,i)+b(13,i)
1018 EE(1,1,-i)= b(10,i)+b(11,i)
1019 EE(2,2,-i)=-b(10,i)+b(11,i)
1020 EE(2,1,-i)=-b(12,i)+b(13,i)
1021 EE(1,2,-i)=-b(12,i)-b(13,i)
1027 c ee(2,1,i)=ee(1,2,i)
1032 write (iout,*) 'Type',i
1034 c write (iout,'(f10.5)') B1(:,i)
1035 write(iout,*) B1(1,i),B1(2,i)
1037 c write (iout,'(f10.5)') B2(:,i)
1038 write(iout,*) B2(1,i),B2(2,i)
1041 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1045 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1049 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1055 C Read electrostatic-interaction parameters
1058 write (iout,'(/a)') 'Electrostatic interaction constants:'
1059 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1060 & 'IT','JT','APP','BPP','AEL6','AEL3'
1062 read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
1063 read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
1064 read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
1065 read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
1070 app (i,j)=epp(i,j)*rri*rri
1071 bpp (i,j)=-2.0D0*epp(i,j)*rri
1072 ael6(i,j)=elpp6(i,j)*4.2D0**6
1073 ael3(i,j)=elpp3(i,j)*4.2D0**3
1075 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1076 & ael6(i,j),ael3(i,j)
1081 C Read side-chain interaction parameters.
1083 read (isidep,*) ipot,expon
1084 if (ipot.lt.1 .or. ipot.gt.5) then
1085 write (iout,'(2a)') 'Error while reading SC interaction',
1086 & 'potential file - unknown potential type.'
1090 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1091 & ', exponents are ',expon,2*expon
1092 goto (10,20,30,30,40) ipot
1093 C----------------------- LJ potential ---------------------------------
1094 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1096 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1097 write (iout,'(a/)') 'The epsilon array:'
1098 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1099 write (iout,'(/a)') 'One-body parameters:'
1100 write (iout,'(a,4x,a)') 'residue','sigma'
1101 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1104 C----------------------- LJK potential --------------------------------
1105 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1106 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1108 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1109 write (iout,'(a/)') 'The epsilon array:'
1110 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1111 write (iout,'(/a)') 'One-body parameters:'
1112 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1113 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1117 C---------------------- GB or BP potential -----------------------------
1119 read (isidep,*)(eps(i,j),j=i,ntyp)
1121 read (isidep,*)(sigma0(i),i=1,ntyp)
1122 read (isidep,*)(sigii(i),i=1,ntyp)
1123 read (isidep,*)(chip(i),i=1,ntyp)
1124 read (isidep,*)(alp(i),i=1,ntyp)
1125 C For the GB potential convert sigma'**2 into chi'
1128 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1132 write (iout,'(/a/)') 'Parameters of the BP potential:'
1133 write (iout,'(a/)') 'The epsilon array:'
1134 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1135 write (iout,'(/a)') 'One-body parameters:'
1136 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1138 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1139 & chip(i),alp(i),i=1,ntyp)
1142 C--------------------- GBV potential -----------------------------------
1143 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1144 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1145 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1147 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1148 write (iout,'(a/)') 'The epsilon array:'
1149 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1150 write (iout,'(/a)') 'One-body parameters:'
1151 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1152 & 's||/s_|_^2',' chip ',' alph '
1153 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1154 & sigii(i),chip(i),alp(i),i=1,ntyp)
1158 C-----------------------------------------------------------------------
1159 C Calculate the "working" parameters of SC interactions.
1167 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1168 sigma(j,i)=sigma(i,j)
1169 rs0(i,j)=dwa16*sigma(i,j)
1173 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1174 & 'Working parameters of the SC interactions:',
1175 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1180 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1189 sigeps=dsign(1.0D0,epsij)
1191 aa(i,j)=epsij*rrij*rrij
1192 bb(i,j)=-sigeps*epsij*rrij
1196 sigt1sq=sigma0(i)**2
1197 sigt2sq=sigma0(j)**2
1200 ratsig1=sigt2sq/sigt1sq
1201 ratsig2=1.0D0/ratsig1
1202 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1203 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1204 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1208 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1209 sigmaii(i,j)=rsum_max
1210 sigmaii(j,i)=rsum_max
1212 c sigmaii(i,j)=r0(i,j)
1213 c sigmaii(j,i)=r0(i,j)
1215 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1216 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1217 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1218 augm(i,j)=epsij*r_augm**(2*expon)
1219 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1226 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1227 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1228 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1233 C Define the SC-p interaction constants
1237 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1239 c aad(i,1)=0.3D0*4.0D0**12
1240 C Following line for constants currently implemented
1241 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1242 aad(i,1)=1.5D0*4.0D0**12
1243 c aad(i,1)=0.17D0*5.6D0**12
1245 C "Soft" SC-p repulsion
1247 C Following line for constants currently implemented
1248 c aad(i,1)=0.3D0*4.0D0**6
1249 C "Hard" SC-p repulsion
1250 bad(i,1)=3.0D0*4.0D0**6
1251 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1260 C 8/9/01 Read the SC-p interaction constants from file
1263 read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1266 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1267 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1268 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1269 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1273 write (iout,*) "Parameters of SC-p interactions:"
1275 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1276 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1281 C Define the constants of the disulfide bridge
1285 c Old arbitrary potential - commented out.
1290 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1291 c energy surface of diethyl disulfide.
1292 c A. Liwo and U. Kozlowska, 11/24/03
1303 write (iout,'(/a)') "Disulfide bridge parameters:"
1304 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1305 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1306 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1307 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,