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)
55 write (iout,*) "iparm",iparm," myparm",myparm
56 c If reading not own parameters, skip assignment
58 if (iparm.eq.myparm .or. .not.separate_parset) then
61 c Setup weights for UNRES
84 call card_concat(controlcard,.false.)
86 c Return if not own parameters
88 if (iparm.ne.myparm .and. separate_parset) return
90 call reads(controlcard,"BONDPAR",bondname_t,bondname)
91 open (ibond,file=bondname_t,status='old')
93 call reads(controlcard,"THETPAR",thetname_t,thetname)
94 open (ithep,file=thetname_t,status='old')
96 call reads(controlcard,"ROTPAR",rotname_t,rotname)
97 open (irotam,file=rotname_t,status='old')
99 call reads(controlcard,"TORPAR",torname_t,torname)
100 open (itorp,file=torname_t,status='old')
102 call reads(controlcard,"TORDPAR",tordname_t,tordname)
103 open (itordp,file=tordname_t,status='old')
105 call reads(controlcard,"SCCORAR",sccorname_t,sccorname)
106 open (isccor,file=sccorname_t,status='old')
108 call reads(controlcard,"FOURIER",fouriername_t,fouriername)
109 open (ifourier,file=fouriername_t,status='old')
111 call reads(controlcard,"ELEPAR",elename_t,elename)
112 open (ielep,file=elename_t,status='old')
114 call reads(controlcard,"SIDEPAR",sidename_t,sidename)
115 open (isidep,file=sidename_t,status='old')
117 call reads(controlcard,"SCPPAR",scpname_t,scpname)
118 open (iscpp,file=scpname_t,status='old')
120 write (iout,*) "Parameter set:",iparm
121 write (iout,*) "Energy-term weights:"
123 write (iout,'(a16,f10.5)') wname(i),ww(i)
125 write (iout,*) "Sidechain potential file : ",
126 & sidename_t(:ilen(sidename_t))
128 write (iout,*) "SCp potential file : ",
129 & scpname_t(:ilen(scpname_t))
131 write (iout,*) "Electrostatic potential file : ",
132 & elename_t(:ilen(elename_t))
133 write (iout,*) "Cumulant coefficient file : ",
134 & fouriername_t(:ilen(fouriername_t))
135 write (iout,*) "Torsional parameter file : ",
136 & torname_t(:ilen(torname_t))
137 write (iout,*) "Double torsional parameter file : ",
138 & tordname_t(:ilen(tordname_t))
139 write (iout,*) "Backbone-rotamer parameter file : ",
140 & sccorname(:ilen(sccorname))
141 write (iout,*) "Bond & inertia constant file : ",
142 & bondname_t(:ilen(bondname_t))
143 write (iout,*) "Bending parameter file : ",
144 & thetname_t(:ilen(thetname_t))
145 write (iout,*) "Rotamer parameter file : ",
146 & rotname_t(:ilen(rotname_t))
149 c Read the virtual-bond parameters, masses, and moments of inertia
150 c and Stokes' radii of the peptide group and side chains
153 read (ibond,*) vbldp0,akp
156 read (ibond,*) vbldsc0(1,i),aksc(1,i)
157 dsc(i) = vbldsc0(1,i)
161 dsc_inv(i)=1.0D0/dsc(i)
165 read (ibond,*) ijunk,vbldp0,akp,rjunk
167 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
169 dsc(i) = vbldsc0(1,i)
173 dsc_inv(i)=1.0D0/dsc(i)
178 write(iout,'(/a/)')"Force constants virtual bonds:"
179 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
181 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
183 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
184 & vbldsc0(1,i),aksc(1,i),abond0(1,i)
186 write (iout,'(13x,3f10.5)')
187 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
193 C Read the parameters of the probability distribution/energy expression
194 C of the virtual-bond valence angles theta
197 read (ithep,*) a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
198 read (ithep,*) (polthet(j,i),j=0,3)
199 read (ithep,*) (gthet(j,i),j=1,3)
200 read (ithep,*) theta0(i),sig0(i),sigc0(i)
206 c & 'Parameters of the virtual-bond valence angles:'
207 c write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
208 c & ' ATHETA0 ',' A1 ',' A2 ',
211 c write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
212 c & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
214 c write (iout,'(/a/9x,5a/79(1h-))')
215 c & 'Parameters of the expression for sigma(theta_c):',
216 c & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
217 c & ' ALPH3 ',' SIGMA0C '
219 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
220 c & (polthet(j,i),j=0,3),sigc0(i)
222 c write (iout,'(/a/9x,5a/79(1h-))')
223 c & 'Parameters of the second gaussian:',
224 c & ' THETA0 ',' SIGMA0 ',' G1 ',
227 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
228 c & sig0(i),(gthet(j,i),j=1,3)
231 & 'Parameters of the virtual-bond valence angles:'
232 write (iout,'(/a/9x,5a/79(1h-))')
233 & 'Coefficients of expansion',
234 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
235 & ' b1*10^1 ',' b2*10^1 '
237 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
238 & a0thet(i),(100*athet(j,i),j=1,2),(10*bthet(j,i),j=1,2)
240 write (iout,'(/a/9x,5a/79(1h-))')
241 & 'Parameters of the expression for sigma(theta_c):',
242 & ' alpha0 ',' alph1 ',' alph2 ',
243 & ' alhp3 ',' sigma0c '
245 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
246 & (polthet(j,i),j=0,3),sigc0(i)
248 write (iout,'(/a/9x,5a/79(1h-))')
249 & 'Parameters of the second gaussian:',
250 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
253 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
254 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
259 C Read the parameters of Utheta determined from ab initio surfaces
260 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
262 read (ithep,*) nthetyp,ntheterm,ntheterm2,
263 & ntheterm3,nsingle,ndouble
264 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
265 read (ithep,*) (ithetyp(i),i=1,ntyp1)
271 aathet(l,i,j,k)=0.0d0
275 bbthet(m,l,i,j,k)=0.0d0
276 ccthet(m,l,i,j,k)=0.0d0
277 ddthet(m,l,i,j,k)=0.0d0
278 eethet(m,l,i,j,k)=0.0d0
284 ffthet(mm,m,l,i,j,k)=0.0d0
285 ggthet(mm,m,l,i,j,k)=0.0d0
295 read (ithep,'(3a)') res1,res2,res3
296 read (ithep,*) aa0thet(i,j,k)
297 read (ithep,*)(aathet(l,i,j,k),l=1,ntheterm)
299 & ((bbthet(lll,ll,i,j,k),lll=1,nsingle),
300 & (ccthet(lll,ll,i,j,k),lll=1,nsingle),
301 & (ddthet(lll,ll,i,j,k),lll=1,nsingle),
302 & (eethet(lll,ll,i,j,k),lll=1,nsingle),ll=1,ntheterm2)
304 & (((ffthet(llll,lll,ll,i,j,k),ffthet(lll,llll,ll,i,j,k),
305 & ggthet(llll,lll,ll,i,j,k),ggthet(lll,llll,ll,i,j,k),
306 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
311 C For dummy ends assign glycine-type coefficients of theta-only terms; the
312 C coefficients of theta-and-gamma-dependent terms are zero.
317 aathet(l,i,j,nthetyp+1)=aathet(l,i,j,1)
318 aathet(l,nthetyp+1,i,j)=aathet(l,1,i,j)
320 aa0thet(i,j,nthetyp+1)=aa0thet(i,j,1)
321 aa0thet(nthetyp+1,i,j)=aa0thet(1,i,j)
324 aathet(l,nthetyp+1,i,nthetyp+1)=aathet(l,1,i,1)
326 aa0thet(nthetyp+1,i,nthetyp+1)=aa0thet(1,i,1)
329 C Control printout of the coefficients of virtual-bond-angle potentials
332 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
336 write (iout,'(//4a)')
337 & 'Type ',onelett(i),onelett(j),onelett(k)
338 write (iout,'(//a,10x,a)') " l","a[l]"
339 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k)
340 write (iout,'(i2,1pe15.5)')
341 & (l,aathet(l,i,j,k),l=1,ntheterm)
343 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
344 & "b",l,"c",l,"d",l,"e",l
346 write (iout,'(i2,4(1pe15.5))') m,
347 & bbthet(m,l,i,j,k),ccthet(m,l,i,j,k),
348 & ddthet(m,l,i,j,k),eethet(m,l,i,j,k)
352 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
353 & "f+",l,"f-",l,"g+",l,"g-",l
356 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
357 & ffthet(n,m,l,i,j,k),ffthet(m,n,l,i,j,k),
358 & ggthet(n,m,l,i,j,k),ggthet(m,n,l,i,j,k)
371 C Read the parameters of the probability distribution/energy expression
372 C of the side chains.
375 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
379 dsc_inv(i)=1.0D0/dsc(i)
390 read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
392 read (irotam,*) bsc(j,i)
393 read (irotam,*) (censc(k,j,i),k=1,3),
394 & ((blower(k,l,j),l=1,k),k=1,3)
401 akl=akl+blower(k,m,j)*blower(l,m,j)
412 write (iout,'(/a)') 'Parameters of side-chain local geometry'
416 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
417 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
418 c write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
419 c write (iout,'(a,f10.4,4(16x,f10.4))')
420 c & 'Center ',(bsc(j,i),j=1,nlobi)
421 c write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
422 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
423 & 'log h',(bsc(j,i),j=1,nlobi)
424 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
425 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
432 c blower(k,l,j)=gaussc(ind,j,i)
437 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
438 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
445 C Read scrot parameters for potentials determined from all-atom AM1 calculations
446 C added by Urszula Kozlowska 07/11/2007
454 read(irotam,*) sc_parmin(j,i)
462 C Read torsional parameters in old format
464 read (itorp,*) ntortyp,nterm_old
465 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
466 read (itorp,*) (itortyp(i),i=1,ntyp)
471 read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
477 write (iout,'(/a/)') 'Torsional constants:'
480 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
481 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
489 C Read torsional parameters
491 read (itorp,*) ntortyp
492 read (itorp,*) (itortyp(i),i=1,ntyp)
493 write (iout,*) 'ntortyp',ntortyp
496 read (itorp,*) nterm(i,j),nlor(i,j)
500 read (itorp,*) kk,v1(k,i,j),v2(k,i,j)
501 v0ij=v0ij+si*v1(k,i,j)
505 read (itorp,*) kk,vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
506 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
513 write (iout,'(/a/)') 'Torsional constants:'
516 write (iout,*) 'ityp',i,' jtyp',j
517 write (iout,*) 'Fourier constants'
519 write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
521 write (iout,*) 'Lorenz constants'
523 write (iout,'(3(1pe15.5))')
524 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
530 C 6/23/01 Read parameters for double torsionals
535 read (itordp,'(3a1)') t1,t2,t3
536 if (t1.ne.onelett(i) .or. t2.ne.onelett(j)
537 & .or. t3.ne.onelett(k)) then
538 write (iout,*) "Error in double torsional parameter file",
540 stop "Error in double torsional parameter file"
542 read (itordp,*) ntermd_1(i,j,k),ntermd_2(i,j,k)
543 read (itordp,*) (v1c(1,l,i,j,k),l=1,ntermd_1(i,j,k))
544 read (itordp,*) (v1s(1,l,i,j,k),l=1,ntermd_1(i,j,k))
545 read (itordp,*) (v1c(2,l,i,j,k),l=1,ntermd_1(i,j,k))
546 read (itordp,*) (v1s(2,l,i,j,k),l=1,ntermd_1(i,j,k))
547 read (itordp,*) ((v2c(l,m,i,j,k),v2c(m,l,i,j,k),
548 & v2s(l,m,i,j,k),v2s(m,l,i,j,k),m=1,l-1),l=1,ntermd_2(i,j,k))
554 write (iout,*) 'Constants for double torsionals'
558 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
559 & ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k)
561 write (iout,*) 'Single angles:'
562 do l=1,ntermd_1(i,j,k)
563 write (iout,'(i5,2f10.5,5x,2f10.5)') l,
564 & v1c(1,l,i,j,k),v1s(1,l,i,j,k),
565 & v1c(2,l,i,j,k),v1s(2,l,i,j,k)
568 write (iout,*) 'Pairs of angles:'
569 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
570 do l=1,ntermd_2(i,j,k)
571 write (iout,'(i5,20f10.5)')
572 & l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k))
575 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
576 do l=1,ntermd_2(i,j,k)
577 write (iout,'(i5,20f10.5)')
578 & l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k))
586 C Read of Side-chain backbone correlation parameters
587 C Modified 11 May 2012 by Adasko
590 read (isccor,*) nsccortyp
591 read (isccor,*) (isccortyp(i),i=1,ntyp)
592 c write (iout,*) 'ntortyp',ntortyp
594 cc maxinter is maximum interaction sites
598 read (isccor,*) nterm_sccor(i,j),nlor_sccor(i,j)
602 do k=1,nterm_sccor(i,j)
603 read (isccor,*) kk,v1sccor(k,l,i,j)
605 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
608 do k=1,nlor_sccor(i,j)
609 read (isccor,*) kk,vlor1sccor(k,i,j),
610 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
611 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
612 &(1+vlor3sccor(k,i,j)**2)
614 v0sccor(i,j)=v0ijsccor
621 write (iout,'(/a/)') 'Torsional constants:'
624 write (iout,*) 'ityp',i,' jtyp',j
625 write (iout,*) 'Fourier constants'
626 do k=1,nterm_sccor(i,j)
627 write(iout,'(2(1pe15.5))')(v1sccor(k,l,i,j),v2sccor(k,l,i,j),
630 write (iout,*) 'Lorenz constants'
631 do k=1,nlor_sccor(i,j)
632 write (iout,'(3(1pe15.5))')
633 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
640 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
641 C interaction energy of the Gly, Ala, and Pro prototypes.
643 read (ifourier,*) nloctyp
646 read (ifourier,*) (b(ii,i),ii=1,13)
648 write (iout,*) 'Type',i
649 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
653 B1tilde(1,i) = b(3,i)
654 B1tilde(2,i) =-b(5,i)
663 Ctilde(2,1,i)=-b(9,i)
671 Dtilde(2,1,i)=-b(8,i)
673 EE(1,1,i)= b(10,i)+b(11,i)
674 EE(2,2,i)=-b(10,i)+b(11,i)
675 EE(2,1,i)= b(12,i)-b(13,i)
676 EE(1,2,i)= b(12,i)+b(13,i)
680 write (iout,*) 'Type',i
682 c write (iout,'(f10.5)') B1(:,i)
683 write(iout,*) B1(1,i),B1(2,i)
685 c write (iout,'(f10.5)') B2(:,i)
686 write(iout,*) B2(1,i),B2(2,i)
689 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
693 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
697 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
702 C Read electrostatic-interaction parameters
705 write (iout,'(/a)') 'Electrostatic interaction constants:'
706 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
707 & 'IT','JT','APP','BPP','AEL6','AEL3'
709 read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
710 read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
711 read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
712 read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
717 app (i,j)=epp(i,j)*rri*rri
718 bpp (i,j)=-2.0D0*epp(i,j)*rri
719 ael6(i,j)=elpp6(i,j)*4.2D0**6
720 ael3(i,j)=elpp3(i,j)*4.2D0**3
721 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
722 & ael6(i,j),ael3(i,j)
726 C Read side-chain interaction parameters.
728 read (isidep,*) ipot,expon
729 if (ipot.lt.1 .or. ipot.gt.6) then
730 write (iout,'(2a)') 'Error while reading SC interaction',
731 & 'potential file - unknown potential type.'
735 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
736 & ', exponents are ',expon,2*expon
737 goto (10,20,30,30,40,50) ipot
738 C----------------------- LJ potential ---------------------------------
739 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
741 write (iout,'(/a/)') 'Parameters of the LJ potential:'
742 write (iout,'(a/)') 'The epsilon array:'
743 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
744 write (iout,'(/a)') 'One-body parameters:'
745 write (iout,'(a,4x,a)') 'residue','sigma'
746 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
749 C----------------------- LJK potential --------------------------------
750 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
751 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
753 write (iout,'(/a/)') 'Parameters of the LJK potential:'
754 write (iout,'(a/)') 'The epsilon array:'
755 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
756 write (iout,'(/a)') 'One-body parameters:'
757 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
758 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
762 C---------------------- GB or BP potential -----------------------------
763 30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
764 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
766 C For the GB potential convert sigma'**2 into chi'
769 chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
773 write (iout,'(/a/)') 'Parameters of the BP potential:'
774 write (iout,'(a/)') 'The epsilon array:'
775 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
776 write (iout,'(/a)') 'One-body parameters:'
777 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
779 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
780 & chip(i),alp(i),i=1,ntyp)
783 C--------------------- GBV potential -----------------------------------
784 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
785 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
786 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
788 write (iout,'(/a/)') 'Parameters of the GBV potential:'
789 write (iout,'(a/)') 'The epsilon array:'
790 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
791 write (iout,'(/a)') 'One-body parameters:'
792 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
793 & 's||/s_|_^2',' chip ',' alph '
794 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
795 & sigii(i),chip(i),alp(i),i=1,ntyp)
798 C--------------------- Momo potential -----------------------------------
802 read (isidep,*) (icharge(i),i=1,ntyp)
803 c write (2,*) "icharge",(icharge(i),i=1,ntyp)
806 c! write (*,*) "Im in ", i, " ", j
808 & eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i),
809 & (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j),
810 & chis(i,j),chis(j,i),
811 & nstate(i,j),(wstate(k,i,j),k=1,4),
816 & dtail(1,i,j),dtail(2,i,j),
817 & epshead(i,j),sig0head(i,j),
818 & rborn(i,j),rborn(j,i),
819 & (wqdip(k,i,j),k=1,2),wquad(i,j),
820 & alphapol(i,j),alphapol(j,i),
821 & (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j)
824 c! write (*,*) "nstate gly-gly", nstate(10,10)
825 c! THIS LOOP FILLS PARAMETERS FOR PAIRS OF AA's NOT EXPLICITLY
826 c! DEFINED IN SCPARM.MOMO. IT DOES SO BY TAKING THEM FROM SYMMETRIC
827 c! PAIR, FOR EG. IF ARG-HIS IS BLANK, IT WILL TAKE PARAMETERS
830 c! DISABLE IT AT >>YOUR OWN PERIL<<
835 sigma(i,j) = sigma(j,i)
836 nstate(i,j) = nstate(j,i)
837 sigmap1(i,j) = sigmap1(j,i)
838 sigmap2(i,j) = sigmap2(j,i)
839 sigiso1(i,j) = sigiso1(j,i)
840 sigiso2(i,j) = sigiso2(j,i)
843 alphasur(k,i,j) = alphasur(k,j,i)
844 wstate(k,i,j) = wstate(k,j,i)
845 alphiso(k,i,j) = alphiso(k,j,i)
848 dhead(2,1,i,j) = dhead(1,1,j,i)
849 dhead(2,2,i,j) = dhead(1,2,j,i)
850 dhead(1,1,i,j) = dhead(2,1,j,i)
851 dhead(1,2,i,j) = dhead(2,2,j,i)
852 dtail(1,i,j) = dtail(1,j,i)
853 dtail(2,i,j) = dtail(2,j,i)
856 c! write (*,*) "dhead(k,l,j,i) = ", dhead(k,l,j,i)
857 c! write (*,*) "dhead(k,l,i,j) = ", dhead(k,l,i,j)
858 c! dhead(l,k,i,j) = dhead(k,l,j,i)
862 epshead(i,j) = epshead(j,i)
863 sig0head(i,j) = sig0head(j,i)
866 wqdip(k,i,j) = wqdip(k,j,i)
869 wquad(i,j) = wquad(j,i)
870 epsintab(i,j) = epsintab(j,i)
875 if (.not.lprint) goto 70
877 & "Parameters of the new physics-based SC-SC interaction potential"
878 write (iout,'(/7a)') 'Residues',' epsGB',' rGB',
879 & ' chi1GB',' chi2GB',' chip1GB',' chip2GB'
882 write (iout,'(2(a3,1x),1pe10.3,5(0pf10.3))')
883 & restyp(i),restyp(j),eps(i,j),sigma(i,j),chi(i,j),chi(j,i),
884 & chipp(i,j),chipp(j,i)
887 write (iout,'(/9a)') 'Residues',' alphasur1',' alphasur2',
888 & ' alphasur3',' alphasur4',' sigmap1',' sigmap2',
892 write (iout,'(2(a3,1x),8(0pf10.3))')
893 & restyp(i),restyp(j),(alphasur(k,i,j),k=1,4),
894 & sigmap1(i,j),sigmap2(j,i),chis(i,j),chis(j,i)
897 write (iout,'(/14a)') 'Residues',' nst ',' wst1',
898 & ' wst2',' wst3',' wst4',' dh11',' dh21',
899 & ' dh12',' dh22',' dt1',' dt2',' epsh1',
903 write (iout,'(2(a3,1x),i3,4f8.3,6f7.2,f9.5,f7.2)')
904 & restyp(i),restyp(j),nstate(i,j),(wstate(k,i,j),k=1,4),
905 & ((dhead(l,k,i,j),l=1,2),k=1,2),dtail(1,i,j),dtail(2,i,j),
906 & epshead(i,j),sig0head(i,j)
909 write (iout,'(/12a)') 'Residues',' ch1',' ch2',
910 & ' rborn1',' rborn2',' wqdip1',' wqdip2',
914 write (iout,'(2(a3,1x),2i4,5f10.3)')
915 & restyp(i),restyp(j),icharge(i),icharge(j),
916 & rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j)
919 write (iout,'(/12a)') 'Residues',
921 & ' alphpol2',' alphiso1',' alpiso2',
922 & ' alpiso3',' alpiso4',' sigiso1',' sigiso2',
926 write (iout,'(2(a3,1x),11f10.3)')
927 & restyp(i),restyp(j),alphapol(i,j),alphapol(j,i),
928 & (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(j,i),
936 C-----------------------------------------------------------------------
937 C Calculate the "working" parameters of SC interactions.
942 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
943 sigma(j,i)=sigma(i,j)
944 rs0(i,j)=dwa16*sigma(i,j)
951 write (iout,*) "IPOT=",ipot
952 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
953 & 'Working parameters of the SC interactions:',
954 & ' a ',' b ',' augm ',' sigma ',' r0 ',
959 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4 .or. ipot.eq.6 ) THEN
968 sigeps=dsign(1.0D0,epsij)
970 aa(i,j)=epsij*rrij*rrij
971 bb(i,j)=-sigeps*epsij*rrij
974 IF ((ipot.gt.2).AND.(ipot.LT.6)) THEN
979 ratsig1=sigt2sq/sigt1sq
980 ratsig2=1.0D0/ratsig1
981 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
982 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
983 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
987 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
988 sigmaii(i,j)=rsum_max
989 sigmaii(j,i)=rsum_max
991 c sigmaii(i,j)=r0(i,j)
992 c sigmaii(j,i)=r0(i,j)
994 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
995 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
996 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
997 augm(i,j)=epsij*r_augm**(2*expon)
998 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1006 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1007 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1008 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1010 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3),2i3,10f8.4,
1012 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1013 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i),
1014 & icharge(i),icharge(j),chipp(i,j),chipp(j,i),
1015 & (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(j,i),
1016 & chis(i,j),chis(j,i),
1017 & nstate(i,j),(wstate(k,i,j),k=1,4),
1018 & ((dhead(l,k,i,j),l=1,2),k=1,2),dtail(1,i,j),dtail(2,i,j),
1019 & epshead(i,j),sig0head(i,j),
1020 & rborn(i,j),(wqdip(k,i,j),k=1,2),wquad(i,j),
1021 & alphapol(i,j),alphapol(j,i),
1022 & (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j)
1030 C Define the SC-p interaction constants
1034 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1036 c aad(i,1)=0.3D0*4.0D0**12
1037 C Following line for constants currently implemented
1038 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1039 aad(i,1)=1.5D0*4.0D0**12
1040 c aad(i,1)=0.17D0*5.6D0**12
1042 C "Soft" SC-p repulsion
1044 C Following line for constants currently implemented
1045 c aad(i,1)=0.3D0*4.0D0**6
1046 C "Hard" SC-p repulsion
1047 bad(i,1)=3.0D0*4.0D0**6
1048 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1057 C 8/9/01 Read the SC-p interaction constants from file
1060 read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1063 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1064 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1065 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1066 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1070 write (iout,*) "Parameters of SC-p interactions:"
1072 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1073 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1078 C Define the constants of the disulfide bridge
1082 c Old arbitrary potential - commented out.
1087 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1088 c energy surface of diethyl disulfide.
1089 c A. Liwo and U. Kozlowska, 11/24/03
1100 write (iout,'(/a)') "Disulfide bridge parameters:"
1101 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1102 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1103 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1104 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,