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 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 cc write (iout,*) "tu dochodze",i
412 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
416 dsc_inv(i)=1.0D0/dsc(i)
427 read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
428 censc(1,1,-i)=censc(1,1,i)
429 censc(2,1,-i)=censc(2,1,i)
430 censc(3,1,-i)=-censc(3,1,i)
432 read (irotam,*) bsc(j,i)
433 read (irotam,*) (censc(k,j,i),k=1,3),
434 & ((blower(k,l,j),l=1,k),k=1,3)
435 censc(1,j,-i)=censc(1,j,i)
436 censc(2,j,-i)=censc(2,j,i)
437 censc(3,j,-i)=-censc(3,j,i)
438 C BSC is amplitude of Gaussian
445 akl=akl+blower(k,m,j)*blower(l,m,j)
449 if (((k.eq.3).and.(l.ne.3))
450 & .or.((l.eq.3).and.(k.ne.3))) then
451 gaussc(k,l,j,-i)=-akl
452 gaussc(l,k,j,-i)=-akl
464 write (iout,'(/a)') 'Parameters of side-chain local geometry'
468 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
469 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
470 c write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
471 c write (iout,'(a,f10.4,4(16x,f10.4))')
472 c & 'Center ',(bsc(j,i),j=1,nlobi)
473 c write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
474 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
475 & 'log h',(bsc(j,i),j=1,nlobi)
476 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
477 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
484 c blower(k,l,j)=gaussc(ind,j,i)
489 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
490 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
497 C Read scrot parameters for potentials determined from all-atom AM1 calculations
498 C added by Urszula Kozlowska 07/11/2007
506 read(irotam,*) sc_parmin(j,i)
514 C Read torsional parameters in old format
516 read (itorp,*) ntortyp,nterm_old
517 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
518 read (itorp,*) (itortyp(i),i=1,ntyp)
523 read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
529 write (iout,'(/a/)') 'Torsional constants:'
532 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
533 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
541 C Read torsional parameters
543 read (itorp,*) ntortyp
544 read (itorp,*) (itortyp(i),i=1,ntyp)
545 write (iout,*) 'ntortyp',ntortyp
548 itortyp(i)=-itortyp(-i)
550 c write (iout,*) 'ntortyp',ntortyp
552 do j=-ntortyp+1,ntortyp-1
553 read (itorp,*) nterm(i,j,iblock),
555 nterm(-i,-j,iblock)=nterm(i,j,iblock)
556 nlor(-i,-j,iblock)=nlor(i,j,iblock)
559 do k=1,nterm(i,j,iblock)
560 read (itorp,*) kk,v1(k,i,j,iblock),
562 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
563 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
564 v0ij=v0ij+si*v1(k,i,j,iblock)
567 do k=1,nlor(i,j,iblock)
568 read (itorp,*) kk,vlor1(k,i,j),
569 & vlor2(k,i,j),vlor3(k,i,j)
570 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
573 v0(-i,-j,iblock)=v0ij
579 write (iout,'(/a/)') 'Torsional constants:'
582 write (iout,*) 'ityp',i,' jtyp',j
583 write (iout,*) 'Fourier constants'
584 do k=1,nterm(i,j,iblock)
585 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
588 write (iout,*) 'Lorenz constants'
589 do k=1,nlor(i,j,iblock)
590 write (iout,'(3(1pe15.5))')
591 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
597 C 6/23/01 Read parameters for double torsionals
601 do j=-ntortyp+1,ntortyp-1
602 do k=-ntortyp+1,ntortyp-1
603 read (itordp,'(3a1)') t1,t2,t3
604 c write (iout,*) "OK onelett",
607 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
608 & .or. t3.ne.toronelet(k)) then
609 write (iout,*) "Error in double torsional parameter file",
612 call MPI_Finalize(Ierror)
614 stop "Error in double torsional parameter file"
616 read (itordp,*) ntermd_1(i,j,k,iblock),
617 & ntermd_2(i,j,k,iblock)
618 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
619 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
620 read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,
621 & ntermd_1(i,j,k,iblock))
622 read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,
623 & ntermd_1(i,j,k,iblock))
624 read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,
625 & ntermd_1(i,j,k,iblock))
626 read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,
627 & ntermd_1(i,j,k,iblock))
628 C Martix of D parameters for one dimesional foureir series
629 do l=1,ntermd_1(i,j,k,iblock)
630 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
631 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
632 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
633 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
634 c write(iout,*) "whcodze" ,
635 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
637 read (itordp,*) ((v2c(l,m,i,j,k,iblock),
638 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
639 & v2s(m,l,i,j,k,iblock),
640 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
641 C Martix of D parameters for two dimesional fourier series
642 do l=1,ntermd_2(i,j,k,iblock)
644 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
645 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
646 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
647 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
656 write (iout,*) 'Constants for double torsionals'
659 do j=-ntortyp+1,ntortyp-1
660 do k=-ntortyp+1,ntortyp-1
661 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
662 & ' nsingle',ntermd_1(i,j,k,iblock),
663 & ' ndouble',ntermd_2(i,j,k,iblock)
665 write (iout,*) 'Single angles:'
666 do l=1,ntermd_1(i,j,k,iblock)
667 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
668 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
669 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
670 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
673 write (iout,*) 'Pairs of angles:'
674 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
675 do l=1,ntermd_2(i,j,k,iblock)
676 write (iout,'(i5,20f10.5)')
677 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
680 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
681 do l=1,ntermd_2(i,j,k,iblock)
682 write (iout,'(i5,20f10.5)')
683 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
684 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
693 C Read of Side-chain backbone correlation parameters
694 C Modified 11 May 2012 by Adasko
697 read (isccor,*) nsccortyp
698 read (isccor,*) (isccortyp(i),i=1,ntyp)
700 isccortyp(i)=-isccortyp(-i)
702 iscprol=isccortyp(20)
703 c write (iout,*) 'ntortyp',ntortyp
705 cc maxinter is maximum interaction sites
710 &nterm_sccor(i,j),nlor_sccor(i,j)
711 write (iout,*) nterm_sccor(i,j)
717 nterm_sccor(-i,j)=nterm_sccor(i,j)
718 nterm_sccor(-i,-j)=nterm_sccor(i,j)
719 nterm_sccor(i,-j)=nterm_sccor(i,j)
720 write (iout,*) nterm_sccor(i,j),nterm_sccor(-i,j),
721 & nterm_sccor(-i,-j),nterm_sccor(i,-j)
722 do k=1,nterm_sccor(i,j)
723 read (isccor,*) kk,v1sccor(k,l,i,j)
725 if (j.eq.iscprol) then
726 if (i.eq.isccortyp(10)) then
727 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
728 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
730 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
731 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
732 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
733 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
734 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
735 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
736 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
737 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
740 if (i.eq.isccortyp(10)) then
741 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
742 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
744 if (j.eq.isccortyp(10)) then
745 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
746 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
748 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
749 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
750 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
751 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
752 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
753 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
757 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
758 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
759 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
760 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
763 do k=1,nlor_sccor(i,j)
764 read (isccor,*) kk,vlor1sccor(k,i,j),
765 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
766 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
767 &(1+vlor3sccor(k,i,j)**2)
769 v0sccor(l,i,j)=v0ijsccor
770 v0sccor(l,-i,j)=v0ijsccor1
771 v0sccor(l,i,-j)=v0ijsccor2
772 v0sccor(l,-i,-j)=v0ijsccor3
778 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
781 write (iout,*) 'ityp',i,' jtyp',j
782 write (iout,*) 'Fourier constants'
783 do k=1,nterm_sccor(i,j)
784 write (iout,'(2(1pe15.5))')
785 & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
787 write (iout,*) 'Lorenz constants'
788 do k=1,nlor_sccor(i,j)
789 write (iout,'(3(1pe15.5))')
790 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
796 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
797 C interaction energy of the Gly, Ala, and Pro prototypes.
799 read (ifourier,*) nloctyp
802 read (ifourier,*) (b(ii,i),ii=1,13)
804 write (iout,*) 'Type',i
805 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
813 B1tilde(1,i) = b(3,i)
814 B1tilde(2,i) =-b(5,i)
815 B1tilde(1,-i) =-b(3,i)
816 B1tilde(2,-i) =b(5,i)
840 Ctilde(2,1,i)=-b(9,i)
842 Ctilde(1,1,-i)=b(7,i)
843 Ctilde(1,2,-i)=-b(9,i)
844 Ctilde(2,1,-i)=b(9,i)
845 Ctilde(2,2,-i)=b(7,i)
847 c Ctilde(1,1,i)=0.0d0
848 c Ctilde(1,2,i)=0.0d0
849 c Ctilde(2,1,i)=0.0d0
850 c Ctilde(2,2,i)=0.0d0
865 Dtilde(2,1,i)=-b(8,i)
867 Dtilde(1,1,-i)=b(6,i)
868 Dtilde(1,2,-i)=-b(8,i)
869 Dtilde(2,1,-i)=b(8,i)
870 Dtilde(2,2,-i)=b(6,i)
872 c Dtilde(1,1,i)=0.0d0
873 c Dtilde(1,2,i)=0.0d0
874 c Dtilde(2,1,i)=0.0d0
875 c Dtilde(2,2,i)=0.0d0
876 EE(1,1,i)= b(10,i)+b(11,i)
877 EE(2,2,i)=-b(10,i)+b(11,i)
878 EE(2,1,i)= b(12,i)-b(13,i)
879 EE(1,2,i)= b(12,i)+b(13,i)
880 EE(1,1,-i)= b(10,i)+b(11,i)
881 EE(2,2,-i)=-b(10,i)+b(11,i)
882 EE(2,1,-i)=-b(12,i)+b(13,i)
883 EE(1,2,-i)=-b(12,i)-b(13,i)
889 c ee(2,1,i)=ee(1,2,i)
894 write (iout,*) 'Type',i
896 c write (iout,'(f10.5)') B1(:,i)
897 write(iout,*) B1(1,i),B1(2,i)
899 c write (iout,'(f10.5)') B2(:,i)
900 write(iout,*) B2(1,i),B2(2,i)
903 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
907 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
911 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
916 C Read electrostatic-interaction parameters
919 write (iout,'(/a)') 'Electrostatic interaction constants:'
920 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
921 & 'IT','JT','APP','BPP','AEL6','AEL3'
923 read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
924 read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
925 read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
926 read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
931 app (i,j)=epp(i,j)*rri*rri
932 bpp (i,j)=-2.0D0*epp(i,j)*rri
933 ael6(i,j)=elpp6(i,j)*4.2D0**6
934 ael3(i,j)=elpp3(i,j)*4.2D0**3
935 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
936 & ael6(i,j),ael3(i,j)
940 C Read side-chain interaction parameters.
942 read (isidep,*) ipot,expon
943 if (ipot.lt.1 .or. ipot.gt.5) then
944 write (iout,'(2a)') 'Error while reading SC interaction',
945 & 'potential file - unknown potential type.'
949 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
950 & ', exponents are ',expon,2*expon
951 goto (10,20,30,30,40) ipot
952 C----------------------- LJ potential ---------------------------------
953 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
955 write (iout,'(/a/)') 'Parameters of the LJ potential:'
956 write (iout,'(a/)') 'The epsilon array:'
957 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
958 write (iout,'(/a)') 'One-body parameters:'
959 write (iout,'(a,4x,a)') 'residue','sigma'
960 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
963 C----------------------- LJK potential --------------------------------
964 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
965 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
967 write (iout,'(/a/)') 'Parameters of the LJK potential:'
968 write (iout,'(a/)') 'The epsilon array:'
969 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
970 write (iout,'(/a)') 'One-body parameters:'
971 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
972 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
976 C---------------------- GB or BP potential -----------------------------
977 30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
978 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
980 C For the GB potential convert sigma'**2 into chi'
983 chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
987 write (iout,'(/a/)') 'Parameters of the BP potential:'
988 write (iout,'(a/)') 'The epsilon array:'
989 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
990 write (iout,'(/a)') 'One-body parameters:'
991 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
993 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
994 & chip(i),alp(i),i=1,ntyp)
997 C--------------------- GBV potential -----------------------------------
998 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
999 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1000 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1002 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1003 write (iout,'(a/)') 'The epsilon array:'
1004 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1005 write (iout,'(/a)') 'One-body parameters:'
1006 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1007 & 's||/s_|_^2',' chip ',' alph '
1008 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1009 & sigii(i),chip(i),alp(i),i=1,ntyp)
1013 C-----------------------------------------------------------------------
1014 C Calculate the "working" parameters of SC interactions.
1022 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1023 sigma(j,i)=sigma(i,j)
1024 rs0(i,j)=dwa16*sigma(i,j)
1028 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1029 & 'Working parameters of the SC interactions:',
1030 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1035 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1044 sigeps=dsign(1.0D0,epsij)
1046 aa(i,j)=epsij*rrij*rrij
1047 bb(i,j)=-sigeps*epsij*rrij
1051 sigt1sq=sigma0(i)**2
1052 sigt2sq=sigma0(j)**2
1055 ratsig1=sigt2sq/sigt1sq
1056 ratsig2=1.0D0/ratsig1
1057 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1058 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1059 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1063 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1064 sigmaii(i,j)=rsum_max
1065 sigmaii(j,i)=rsum_max
1067 c sigmaii(i,j)=r0(i,j)
1068 c sigmaii(j,i)=r0(i,j)
1070 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1071 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1072 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1073 augm(i,j)=epsij*r_augm**(2*expon)
1074 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1081 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1082 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1083 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1088 C Define the SC-p interaction constants
1092 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1094 c aad(i,1)=0.3D0*4.0D0**12
1095 C Following line for constants currently implemented
1096 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1097 aad(i,1)=1.5D0*4.0D0**12
1098 c aad(i,1)=0.17D0*5.6D0**12
1100 C "Soft" SC-p repulsion
1102 C Following line for constants currently implemented
1103 c aad(i,1)=0.3D0*4.0D0**6
1104 C "Hard" SC-p repulsion
1105 bad(i,1)=3.0D0*4.0D0**6
1106 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1115 C 8/9/01 Read the SC-p interaction constants from file
1118 read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1121 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1122 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1123 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1124 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1128 write (iout,*) "Parameters of SC-p interactions:"
1130 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1131 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1136 C Define the constants of the disulfide bridge
1140 c Old arbitrary potential - commented out.
1145 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1146 c energy surface of diethyl disulfide.
1147 c A. Liwo and U. Kozlowska, 11/24/03
1158 write (iout,'(/a)') "Disulfide bridge parameters:"
1159 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1160 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1161 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1162 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,