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
39 C write (iout,*) "KURWA"
43 C Set LPRINT=.TRUE. for debugging
44 dwa16=2.0d0**(1.0d0/6.0d0)
47 C Assign virtual-bond length
51 call card_concat(controlcard,.true.)
54 key = wname(i)(:ilen(wname(i)))
55 call reada(controlcard,key(:ilen(key)),ww(i),1.0d0)
58 write (iout,*) "iparm",iparm," myparm",myparm
59 c If reading not own parameters, skip assignment
60 call reada(controlcard,"D0CM",d0cm,3.78d0)
61 call reada(controlcard,"AKCM",akcm,15.1d0)
62 call reada(controlcard,"AKTH",akth,11.0d0)
63 call reada(controlcard,"AKCT",akct,12.0d0)
64 call reada(controlcard,"V1SS",v1ss,-1.08d0)
65 call reada(controlcard,"V2SS",v2ss,7.61d0)
66 call reada(controlcard,"V3SS",v3ss,13.7d0)
67 call reada(controlcard,"EBR",ebr,-5.50D0)
68 call reada(controlcard,"DTRISS",dtriss,1.0D0)
69 call reada(controlcard,"ATRISS",atriss,0.3D0)
70 call reada(controlcard,"BTRISS",btriss,0.02D0)
71 call reada(controlcard,"CTRISS",ctriss,1.0D0)
72 dyn_ss=(index(controlcard,'DYN_SS').gt.0)
74 C dyn_ss_mask(i)=.false.
78 c Old arbitrary potential - commented out.
83 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
84 c energy surface of diethyl disulfide.
85 c A. Liwo and U. Kozlowska, 11/24/03
97 dyn_ssbond_ij(i,j)=1.0d300
100 call reada(controlcard,"HT",Ht,0.0D0)
102 C ss_depth=ebr/wsc-0.25*eps(1,1)
103 C write(iout,*) HT,wsc,eps(1,1),'KURWA'
104 C Ht=Ht/wsc-0.25*eps(1,1)
113 C ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
116 if (iparm.eq.myparm .or. .not.separate_parset) then
119 c Setup weights for UNRES
141 call card_concat(controlcard,.false.)
143 c Return if not own parameters
145 if (iparm.ne.myparm .and. separate_parset) return
147 call reads(controlcard,"BONDPAR",bondname_t,bondname)
148 open (ibond,file=bondname_t,status='old')
150 call reads(controlcard,"THETPAR",thetname_t,thetname)
151 open (ithep,file=thetname_t,status='old')
153 call reads(controlcard,"ROTPAR",rotname_t,rotname)
154 open (irotam,file=rotname_t,status='old')
156 call reads(controlcard,"TORPAR",torname_t,torname)
157 open (itorp,file=torname_t,status='old')
159 call reads(controlcard,"TORDPAR",tordname_t,tordname)
160 open (itordp,file=tordname_t,status='old')
162 call reads(controlcard,"SCCORPAR",sccorname_t,sccorname)
163 open (isccor,file=sccorname_t,status='old')
165 call reads(controlcard,"FOURIER",fouriername_t,fouriername)
166 open (ifourier,file=fouriername_t,status='old')
168 call reads(controlcard,"ELEPAR",elename_t,elename)
169 open (ielep,file=elename_t,status='old')
171 call reads(controlcard,"SIDEPAR",sidename_t,sidename)
172 open (isidep,file=sidename_t,status='old')
174 call reads(controlcard,"SCPPAR",scpname_t,scpname)
175 open (iscpp,file=scpname_t,status='old')
177 write (iout,*) "Parameter set:",iparm
178 write (iout,*) "Energy-term weights:"
180 write (iout,'(a16,f10.5)') wname(i),ww(i)
182 write (iout,*) "Sidechain potential file : ",
183 & sidename_t(:ilen(sidename_t))
185 write (iout,*) "SCp potential file : ",
186 & scpname_t(:ilen(scpname_t))
188 write (iout,*) "Electrostatic potential file : ",
189 & elename_t(:ilen(elename_t))
190 write (iout,*) "Cumulant coefficient file : ",
191 & fouriername_t(:ilen(fouriername_t))
192 write (iout,*) "Torsional parameter file : ",
193 & torname_t(:ilen(torname_t))
194 write (iout,*) "Double torsional parameter file : ",
195 & tordname_t(:ilen(tordname_t))
196 write (iout,*) "Backbone-rotamer parameter file : ",
197 & sccorname(:ilen(sccorname))
198 write (iout,*) "Bond & inertia constant file : ",
199 & bondname_t(:ilen(bondname_t))
200 write (iout,*) "Bending parameter file : ",
201 & thetname_t(:ilen(thetname_t))
202 write (iout,*) "Rotamer parameter file : ",
203 & rotname_t(:ilen(rotname_t))
206 c Read the virtual-bond parameters, masses, and moments of inertia
207 c and Stokes' radii of the peptide group and side chains
210 read (ibond,*) vbldp0,akp
213 read (ibond,*) vbldsc0(1,i),aksc(1,i)
214 dsc(i) = vbldsc0(1,i)
218 dsc_inv(i)=1.0D0/dsc(i)
222 read (ibond,*) ijunk,vbldp0,akp,rjunk
224 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
226 dsc(i) = vbldsc0(1,i)
230 dsc_inv(i)=1.0D0/dsc(i)
235 write(iout,'(/a/)')"Force constants virtual bonds:"
236 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
238 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
240 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
241 & vbldsc0(1,i),aksc(1,i),abond0(1,i)
243 write (iout,'(13x,3f10.5)')
244 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
250 C Read the parameters of the probability distribution/energy expression
251 C of the virtual-bond valence angles theta
254 read (ithep,*) a0thet(i),(athet(j,i,1,1),j=1,2),
255 & (bthet(j,i,1,1),j=1,2)
256 read (ithep,*) (polthet(j,i),j=0,3)
257 read (ithep,*) (gthet(j,i),j=1,3)
258 read (ithep,*) theta0(i),sig0(i),sigc0(i)
262 athet(1,i,1,-1)=athet(1,i,1,1)
263 athet(2,i,1,-1)=athet(2,i,1,1)
264 bthet(1,i,1,-1)=-bthet(1,i,1,1)
265 bthet(2,i,1,-1)=-bthet(2,i,1,1)
266 athet(1,i,-1,1)=-athet(1,i,1,1)
267 athet(2,i,-1,1)=-athet(2,i,1,1)
268 bthet(1,i,-1,1)=bthet(1,i,1,1)
269 bthet(2,i,-1,1)=bthet(2,i,1,1)
273 athet(1,i,-1,-1)=athet(1,-i,1,1)
274 athet(2,i,-1,-1)=-athet(2,-i,1,1)
275 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
276 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
277 athet(1,i,-1,1)=athet(1,-i,1,1)
278 athet(2,i,-1,1)=-athet(2,-i,1,1)
279 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
280 bthet(2,i,-1,1)=bthet(2,-i,1,1)
281 athet(1,i,1,-1)=-athet(1,-i,1,1)
282 athet(2,i,1,-1)=athet(2,-i,1,1)
283 bthet(1,i,1,-1)=bthet(1,-i,1,1)
284 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
289 polthet(j,i)=polthet(j,-i)
292 gthet(j,i)=gthet(j,-i)
298 c & 'Parameters of the virtual-bond valence angles:'
299 c write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
300 c & ' ATHETA0 ',' A1 ',' A2 ',
303 c write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
304 c & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
306 c write (iout,'(/a/9x,5a/79(1h-))')
307 c & 'Parameters of the expression for sigma(theta_c):',
308 c & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
309 c & ' ALPH3 ',' SIGMA0C '
311 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
312 c & (polthet(j,i),j=0,3),sigc0(i)
314 c write (iout,'(/a/9x,5a/79(1h-))')
315 c & 'Parameters of the second gaussian:',
316 c & ' THETA0 ',' SIGMA0 ',' G1 ',
319 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
320 c & sig0(i),(gthet(j,i),j=1,3)
323 & 'Parameters of the virtual-bond valence angles:'
324 write (iout,'(/a/9x,5a/79(1h-))')
325 & 'Coefficients of expansion',
326 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
327 & ' b1*10^1 ',' b2*10^1 '
329 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
330 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
331 & (10*bthet(j,i,1,1),j=1,2)
333 write (iout,'(/a/9x,5a/79(1h-))')
334 & 'Parameters of the expression for sigma(theta_c):',
335 & ' alpha0 ',' alph1 ',' alph2 ',
336 & ' alhp3 ',' sigma0c '
338 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
339 & (polthet(j,i),j=0,3),sigc0(i)
341 write (iout,'(/a/9x,5a/79(1h-))')
342 & 'Parameters of the second gaussian:',
343 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
346 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
347 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
352 C Read the parameters of Utheta determined from ab initio surfaces
353 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
355 write (iout,*) "tu dochodze"
356 read (ithep,*) nthetyp,ntheterm,ntheterm2,
357 & ntheterm3,nsingle,ndouble
358 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
359 read (ithep,*) (ithetyp(i),i=1,ntyp1)
361 ithetyp(i)=-ithetyp(-i)
363 write (iout,*) "tu dochodze"
365 do i=-maxthetyp,maxthetyp
366 do j=-maxthetyp,maxthetyp
367 do k=-maxthetyp,maxthetyp
368 aa0thet(i,j,k,iblock)=0.0d0
370 aathet(l,i,j,k,iblock)=0.0d0
374 bbthet(m,l,i,j,k,iblock)=0.0d0
375 ccthet(m,l,i,j,k,iblock)=0.0d0
376 ddthet(m,l,i,j,k,iblock)=0.0d0
377 eethet(m,l,i,j,k,iblock)=0.0d0
383 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
384 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
392 write (iout,*) "KURWA1"
395 do j=-nthetyp,nthetyp
396 do k=-nthetyp,nthetyp
397 read (ithep,'(6a)') res1
398 write(iout,*) res1,i,j,k
399 read (ithep,*) aa0thet(i,j,k,iblock)
400 read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
402 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
403 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
404 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
405 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle)
408 & (((ffthet(llll,lll,ll,i,j,k,iblock),
409 & ffthet(lll,llll,ll,i,j,k,iblock),
410 & ggthet(llll,lll,ll,i,j,k,iblock)
411 & ,ggthet(lll,llll,ll,i,j,k,iblock),
412 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
416 write(iout,*) "KURWA1.1"
418 C For dummy ends assign glycine-type coefficients of theta-only terms; the
419 C coefficients of theta-and-gamma-dependent terms are zero.
424 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
425 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
427 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
428 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
431 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
433 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
436 write(iout,*) "KURWA1.5"
437 C Substitution for D aminoacids from symmetry.
440 do j=-nthetyp,nthetyp
441 do k=-nthetyp,nthetyp
442 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
444 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
448 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
449 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
450 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
451 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
457 ffthet(llll,lll,ll,i,j,k,iblock)=
458 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
459 ffthet(lll,llll,ll,i,j,k,iblock)=
460 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
461 ggthet(llll,lll,ll,i,j,k,iblock)=
462 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
463 ggthet(lll,llll,ll,i,j,k,iblock)=
464 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
474 C Control printout of the coefficients of virtual-bond-angle potentials
477 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
481 write (iout,'(//4a)')
482 & 'Type ',onelett(i),onelett(j),onelett(k)
483 write (iout,'(//a,10x,a)') " l","a[l]"
484 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
485 write (iout,'(i2,1pe15.5)')
486 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
488 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
489 & "b",l,"c",l,"d",l,"e",l
491 write (iout,'(i2,4(1pe15.5))') m,
492 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
493 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
497 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
498 & "f+",l,"f-",l,"g+",l,"g-",l
501 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
502 & ffthet(n,m,l,i,j,k,iblock),
503 & ffthet(m,n,l,i,j,k,iblock),
504 & ggthet(n,m,l,i,j,k,iblock),
505 & ggthet(m,n,l,i,j,k,iblock)
515 write(iout,*) 'KURWA2'
518 C Read the parameters of the probability distribution/energy expression
519 C of the side chains.
522 cc write (iout,*) "tu dochodze",i
523 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
527 dsc_inv(i)=1.0D0/dsc(i)
538 read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
539 censc(1,1,-i)=censc(1,1,i)
540 censc(2,1,-i)=censc(2,1,i)
541 censc(3,1,-i)=-censc(3,1,i)
543 read (irotam,*) bsc(j,i)
544 read (irotam,*) (censc(k,j,i),k=1,3),
545 & ((blower(k,l,j),l=1,k),k=1,3)
546 censc(1,j,-i)=censc(1,j,i)
547 censc(2,j,-i)=censc(2,j,i)
548 censc(3,j,-i)=-censc(3,j,i)
549 C BSC is amplitude of Gaussian
556 akl=akl+blower(k,m,j)*blower(l,m,j)
560 if (((k.eq.3).and.(l.ne.3))
561 & .or.((l.eq.3).and.(k.ne.3))) then
562 gaussc(k,l,j,-i)=-akl
563 gaussc(l,k,j,-i)=-akl
575 write (iout,'(/a)') 'Parameters of side-chain local geometry'
579 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
580 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
581 c write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
582 c write (iout,'(a,f10.4,4(16x,f10.4))')
583 c & 'Center ',(bsc(j,i),j=1,nlobi)
584 c write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
585 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
586 & 'log h',(bsc(j,i),j=1,nlobi)
587 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
588 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
595 c blower(k,l,j)=gaussc(ind,j,i)
600 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
601 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
608 C Read scrot parameters for potentials determined from all-atom AM1 calculations
609 C added by Urszula Kozlowska 07/11/2007
617 read(irotam,*) sc_parmin(j,i)
623 write (iout,*) 'KURWAKURWA'
626 C Read torsional parameters in old format
628 read (itorp,*) ntortyp,nterm_old
629 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
630 read (itorp,*) (itortyp(i),i=1,ntyp)
635 read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
641 write (iout,'(/a/)') 'Torsional constants:'
644 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
645 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
653 C Read torsional parameters
655 read (itorp,*) ntortyp
656 read (itorp,*) (itortyp(i),i=1,ntyp)
657 write (iout,*) 'ntortyp',ntortyp
660 itortyp(i)=-itortyp(-i)
662 c write (iout,*) 'ntortyp',ntortyp
664 do j=-ntortyp+1,ntortyp-1
665 read (itorp,*) nterm(i,j,iblock),
667 nterm(-i,-j,iblock)=nterm(i,j,iblock)
668 nlor(-i,-j,iblock)=nlor(i,j,iblock)
671 do k=1,nterm(i,j,iblock)
672 read (itorp,*) kk,v1(k,i,j,iblock),
674 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
675 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
676 v0ij=v0ij+si*v1(k,i,j,iblock)
679 do k=1,nlor(i,j,iblock)
680 read (itorp,*) kk,vlor1(k,i,j),
681 & vlor2(k,i,j),vlor3(k,i,j)
682 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
685 v0(-i,-j,iblock)=v0ij
691 write (iout,'(/a/)') 'Torsional constants:'
694 write (iout,*) 'ityp',i,' jtyp',j
695 write (iout,*) 'Fourier constants'
696 do k=1,nterm(i,j,iblock)
697 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
700 write (iout,*) 'Lorenz constants'
701 do k=1,nlor(i,j,iblock)
702 write (iout,'(3(1pe15.5))')
703 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
709 C 6/23/01 Read parameters for double torsionals
713 do j=-ntortyp+1,ntortyp-1
714 do k=-ntortyp+1,ntortyp-1
715 read (itordp,'(3a1)') t1,t2,t3
716 c write (iout,*) "OK onelett",
719 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
720 & .or. t3.ne.toronelet(k)) then
721 write (iout,*) "Error in double torsional parameter file",
724 call MPI_Finalize(Ierror)
726 stop "Error in double torsional parameter file"
728 read (itordp,*) ntermd_1(i,j,k,iblock),
729 & ntermd_2(i,j,k,iblock)
730 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
731 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
732 read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,
733 & ntermd_1(i,j,k,iblock))
734 read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,
735 & ntermd_1(i,j,k,iblock))
736 read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,
737 & ntermd_1(i,j,k,iblock))
738 read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,
739 & ntermd_1(i,j,k,iblock))
740 C Martix of D parameters for one dimesional foureir series
741 do l=1,ntermd_1(i,j,k,iblock)
742 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
743 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
744 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
745 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
746 c write(iout,*) "whcodze" ,
747 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
749 read (itordp,*) ((v2c(l,m,i,j,k,iblock),
750 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
751 & v2s(m,l,i,j,k,iblock),
752 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
753 C Martix of D parameters for two dimesional fourier series
754 do l=1,ntermd_2(i,j,k,iblock)
756 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
757 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
758 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
759 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
768 write (iout,*) 'Constants for double torsionals'
771 do j=-ntortyp+1,ntortyp-1
772 do k=-ntortyp+1,ntortyp-1
773 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
774 & ' nsingle',ntermd_1(i,j,k,iblock),
775 & ' ndouble',ntermd_2(i,j,k,iblock)
777 write (iout,*) 'Single angles:'
778 do l=1,ntermd_1(i,j,k,iblock)
779 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
780 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
781 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
782 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
785 write (iout,*) 'Pairs of angles:'
786 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
787 do l=1,ntermd_2(i,j,k,iblock)
788 write (iout,'(i5,20f10.5)')
789 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
792 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
793 do l=1,ntermd_2(i,j,k,iblock)
794 write (iout,'(i5,20f10.5)')
795 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
796 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
805 C Read of Side-chain backbone correlation parameters
806 C Modified 11 May 2012 by Adasko
809 read (isccor,*) nsccortyp
810 read (isccor,*) (isccortyp(i),i=1,ntyp)
812 isccortyp(i)=-isccortyp(-i)
814 iscprol=isccortyp(20)
815 c write (iout,*) 'ntortyp',ntortyp
817 cc maxinter is maximum interaction sites
822 &nterm_sccor(i,j),nlor_sccor(i,j)
823 write (iout,*) nterm_sccor(i,j)
829 nterm_sccor(-i,j)=nterm_sccor(i,j)
830 nterm_sccor(-i,-j)=nterm_sccor(i,j)
831 nterm_sccor(i,-j)=nterm_sccor(i,j)
832 write (iout,*) nterm_sccor(i,j),nterm_sccor(-i,j),
833 & nterm_sccor(-i,-j),nterm_sccor(i,-j)
834 do k=1,nterm_sccor(i,j)
835 read (isccor,*) kk,v1sccor(k,l,i,j)
837 if (j.eq.iscprol) then
838 if (i.eq.isccortyp(10)) then
839 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
840 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
842 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
843 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
844 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
845 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
846 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
847 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
848 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
849 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
852 if (i.eq.isccortyp(10)) then
853 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
854 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
856 if (j.eq.isccortyp(10)) then
857 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
858 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
860 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
861 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
862 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
863 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
864 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
865 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
869 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
870 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
871 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
872 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
875 do k=1,nlor_sccor(i,j)
876 read (isccor,*) kk,vlor1sccor(k,i,j),
877 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
878 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
879 &(1+vlor3sccor(k,i,j)**2)
881 v0sccor(l,i,j)=v0ijsccor
882 v0sccor(l,-i,j)=v0ijsccor1
883 v0sccor(l,i,-j)=v0ijsccor2
884 v0sccor(l,-i,-j)=v0ijsccor3
890 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
893 write (iout,*) 'ityp',i,' jtyp',j
894 write (iout,*) 'Fourier constants'
895 do k=1,nterm_sccor(i,j)
896 write (iout,'(2(1pe15.5))')
897 & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
899 write (iout,*) 'Lorenz constants'
900 do k=1,nlor_sccor(i,j)
901 write (iout,'(3(1pe15.5))')
902 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
908 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
909 C interaction energy of the Gly, Ala, and Pro prototypes.
911 read (ifourier,*) nloctyp
914 read (ifourier,*) (b(ii,i),ii=1,13)
916 write (iout,*) 'Type',i
917 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
925 B1tilde(1,i) = b(3,i)
926 B1tilde(2,i) =-b(5,i)
927 B1tilde(1,-i) =-b(3,i)
928 B1tilde(2,-i) =b(5,i)
952 Ctilde(2,1,i)=-b(9,i)
954 Ctilde(1,1,-i)=b(7,i)
955 Ctilde(1,2,-i)=-b(9,i)
956 Ctilde(2,1,-i)=b(9,i)
957 Ctilde(2,2,-i)=b(7,i)
959 c Ctilde(1,1,i)=0.0d0
960 c Ctilde(1,2,i)=0.0d0
961 c Ctilde(2,1,i)=0.0d0
962 c Ctilde(2,2,i)=0.0d0
977 Dtilde(2,1,i)=-b(8,i)
979 Dtilde(1,1,-i)=b(6,i)
980 Dtilde(1,2,-i)=-b(8,i)
981 Dtilde(2,1,-i)=b(8,i)
982 Dtilde(2,2,-i)=b(6,i)
984 c Dtilde(1,1,i)=0.0d0
985 c Dtilde(1,2,i)=0.0d0
986 c Dtilde(2,1,i)=0.0d0
987 c Dtilde(2,2,i)=0.0d0
988 EE(1,1,i)= b(10,i)+b(11,i)
989 EE(2,2,i)=-b(10,i)+b(11,i)
990 EE(2,1,i)= b(12,i)-b(13,i)
991 EE(1,2,i)= b(12,i)+b(13,i)
992 EE(1,1,-i)= b(10,i)+b(11,i)
993 EE(2,2,-i)=-b(10,i)+b(11,i)
994 EE(2,1,-i)=-b(12,i)+b(13,i)
995 EE(1,2,-i)=-b(12,i)-b(13,i)
1001 c ee(2,1,i)=ee(1,2,i)
1006 write (iout,*) 'Type',i
1008 c write (iout,'(f10.5)') B1(:,i)
1009 write(iout,*) B1(1,i),B1(2,i)
1011 c write (iout,'(f10.5)') B2(:,i)
1012 write(iout,*) B2(1,i),B2(2,i)
1015 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1019 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1023 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1028 C Read electrostatic-interaction parameters
1031 write (iout,'(/a)') 'Electrostatic interaction constants:'
1032 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1033 & 'IT','JT','APP','BPP','AEL6','AEL3'
1035 read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
1036 read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
1037 read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
1038 read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
1043 app (i,j)=epp(i,j)*rri*rri
1044 bpp (i,j)=-2.0D0*epp(i,j)*rri
1045 ael6(i,j)=elpp6(i,j)*4.2D0**6
1046 ael3(i,j)=elpp3(i,j)*4.2D0**3
1048 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1049 & ael6(i,j),ael3(i,j)
1054 C Read side-chain interaction parameters.
1056 read (isidep,*) ipot,expon
1057 if (ipot.lt.1 .or. ipot.gt.5) then
1058 write (iout,'(2a)') 'Error while reading SC interaction',
1059 & 'potential file - unknown potential type.'
1063 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1064 & ', exponents are ',expon,2*expon
1065 goto (10,20,30,30,40) ipot
1066 C----------------------- LJ potential ---------------------------------
1067 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1069 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1070 write (iout,'(a/)') 'The epsilon array:'
1071 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1072 write (iout,'(/a)') 'One-body parameters:'
1073 write (iout,'(a,4x,a)') 'residue','sigma'
1074 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1077 C----------------------- LJK potential --------------------------------
1078 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1079 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1081 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1082 write (iout,'(a/)') 'The epsilon array:'
1083 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1084 write (iout,'(/a)') 'One-body parameters:'
1085 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1086 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1090 C---------------------- GB or BP potential -----------------------------
1092 read (isidep,*)(eps(i,j),j=i,ntyp)
1094 read (isidep,*)(sigma0(i),i=1,ntyp)
1095 read (isidep,*)(sigii(i),i=1,ntyp)
1096 read (isidep,*)(chip(i),i=1,ntyp)
1097 read (isidep,*)(alp(i),i=1,ntyp)
1098 C For the GB potential convert sigma'**2 into chi'
1101 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1105 write (iout,'(/a/)') 'Parameters of the BP potential:'
1106 write (iout,'(a/)') 'The epsilon array:'
1107 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1108 write (iout,'(/a)') 'One-body parameters:'
1109 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1111 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1112 & chip(i),alp(i),i=1,ntyp)
1115 C--------------------- GBV potential -----------------------------------
1116 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1117 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1118 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1120 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1121 write (iout,'(a/)') 'The epsilon array:'
1122 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1123 write (iout,'(/a)') 'One-body parameters:'
1124 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1125 & 's||/s_|_^2',' chip ',' alph '
1126 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1127 & sigii(i),chip(i),alp(i),i=1,ntyp)
1131 C-----------------------------------------------------------------------
1132 C Calculate the "working" parameters of SC interactions.
1140 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1141 sigma(j,i)=sigma(i,j)
1142 rs0(i,j)=dwa16*sigma(i,j)
1146 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1147 & 'Working parameters of the SC interactions:',
1148 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1153 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1162 sigeps=dsign(1.0D0,epsij)
1164 aa(i,j)=epsij*rrij*rrij
1165 bb(i,j)=-sigeps*epsij*rrij
1169 sigt1sq=sigma0(i)**2
1170 sigt2sq=sigma0(j)**2
1173 ratsig1=sigt2sq/sigt1sq
1174 ratsig2=1.0D0/ratsig1
1175 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1176 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1177 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1181 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1182 sigmaii(i,j)=rsum_max
1183 sigmaii(j,i)=rsum_max
1185 c sigmaii(i,j)=r0(i,j)
1186 c sigmaii(j,i)=r0(i,j)
1188 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1189 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1190 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1191 augm(i,j)=epsij*r_augm**(2*expon)
1192 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1199 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1200 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1201 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1206 C Define the SC-p interaction constants
1210 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1212 c aad(i,1)=0.3D0*4.0D0**12
1213 C Following line for constants currently implemented
1214 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1215 aad(i,1)=1.5D0*4.0D0**12
1216 c aad(i,1)=0.17D0*5.6D0**12
1218 C "Soft" SC-p repulsion
1220 C Following line for constants currently implemented
1221 c aad(i,1)=0.3D0*4.0D0**6
1222 C "Hard" SC-p repulsion
1223 bad(i,1)=3.0D0*4.0D0**6
1224 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1233 C 8/9/01 Read the SC-p interaction constants from file
1236 read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1239 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1240 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1241 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1242 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1246 write (iout,*) "Parameters of SC-p interactions:"
1248 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1249 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1254 C Define the constants of the disulfide bridge
1258 c Old arbitrary potential - commented out.
1263 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1264 c energy surface of diethyl disulfide.
1265 c A. Liwo and U. Kozlowska, 11/24/03
1274 write (iout,*) dyn_ss,'dyndyn'
1276 ss_depth=ebr/wsc-0.25*eps(1,1)
1277 write(iout,*) akcm,whpb,wsc,'KURWA'
1278 Ht=Ht/wsc-0.25*eps(1,1)
1287 ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
1291 write (iout,'(/a)') "Disulfide bridge parameters:"
1292 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1293 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1294 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1295 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,