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 'COMMON.IOUNITS'
10 include 'COMMON.CHAIN'
11 include 'COMMON.INTERACT'
13 include 'COMMON.LOCAL'
14 include 'COMMON.TORSION'
15 include 'COMMON.FFIELD'
16 include 'COMMON.NAMES'
17 include 'COMMON.SBRIDGE'
18 include 'COMMON.WEIGHTS'
19 include 'COMMON.SCCOR'
20 include 'COMMON.SCROT'
21 include 'COMMON.SHIELD'
22 include 'COMMON.CONTROL'
23 include 'COMMON.DERIV'
24 include 'COMMON.OPTIM'
25 include 'COMMON.ENERGIES'
27 character*1 onelett(4) /"G","A","P","D"/
28 character*1 toronelet(-2:2) /"p","a","G","A","P"/
30 dimension blower(3,3,maxlob)
31 character*800 controlcard
32 character*256 bondname_t,thetname_t,rotname_t,torname_t,
33 & tordname_t,fouriername_t,elename_t,sidename_t,scpname_t,
39 double precision ip,mp
41 character*3 lancuch,ucase
45 call getenv("PRINT_PARM",lancuch)
46 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
47 write (iout,*) "lprint ",lprint
48 C Set LPRINT=.TRUE. for debugging
49 dwa16=2.0d0**(1.0d0/6.0d0)
51 C Assign virtual-bond length
55 call card_concat(controlcard,.true.)
58 key = wname(i)(:ilen(wname(i)))
59 call reada(controlcard,key(:ilen(key)),ww(i),0.0d0)
61 c print *,me," checkpoint 7"
63 write (iout,*) "iparm",iparm," myparm",myparm
65 c If reading not own parameters, skip assignment
66 call reada(controlcard,"D0CM",d0cm,3.78d0)
67 call reada(controlcard,"AKCM",akcm,15.1d0)
68 call reada(controlcard,"AKTH",akth,11.0d0)
69 call reada(controlcard,"AKCT",akct,12.0d0)
70 call reada(controlcard,"V1SS",v1ss,-1.08d0)
71 call reada(controlcard,"V2SS",v2ss,7.61d0)
72 call reada(controlcard,"V3SS",v3ss,13.7d0)
73 call reada(controlcard,"EBR",ebr,-5.50D0)
74 call reada(controlcard,"DTRISS",dtriss,1.0D0)
75 call reada(controlcard,"ATRISS",atriss,0.3D0)
76 call reada(controlcard,"BTRISS",btriss,0.02D0)
77 call reada(controlcard,"CTRISS",ctriss,1.0D0)
78 dyn_ss=(index(controlcard,'DYN_SS').gt.0)
79 write(iout,*) "ATRISS",atriss
80 write(iout,*) "BTRISS",btriss
81 write(iout,*) "CTRISS",ctriss
82 write(iout,*) "DTRISS",dtriss
85 C dyn_ss_mask(i)=.false.
89 c Old arbitrary potential - commented out.
94 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
95 c energy surface of diethyl disulfide.
96 c A. Liwo and U. Kozlowska, 11/24/03
108 dyn_ssbond_ij(i,j)=1.0d300
111 call reada(controlcard,"HT",Ht,0.0D0)
113 C ss_depth=ebr/wsc-0.25*eps(1,1)
114 C write(iout,*) HT,wsc,eps(1,1),'KURWA'
115 C Ht=Ht/wsc-0.25*eps(1,1)
124 C ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
127 if (iparm.eq.myparm) then
130 c Setup weights for UNRES
155 call card_concat(controlcard,.false.)
157 c Return if not own parameters
159 if (iparm.ne.myparm) return
161 call reads(controlcard,"BONDPAR",bondname_t,bondname)
162 open (ibond,file=bondname_t,status='old')
164 call reads(controlcard,"THETPAR",thetname_t,thetname)
165 open (ithep,file=thetname_t,status='old')
167 call reads(controlcard,"ROTPAR",rotname_t,rotname)
168 open (irotam,file=rotname_t,status='old')
170 call reads(controlcard,"TORPAR",torname_t,torname)
171 open (itorp,file=torname_t,status='old')
173 call reads(controlcard,"TORDPAR",tordname_t,tordname)
174 write (iout,*) "tor_mode",tor_mode
177 & open (itordp,file=tordname_t,status='old')
179 call reads(controlcard,"SCCORPAR",sccorname_t,sccorname)
180 open (isccor,file=sccorname_t,status='old')
182 call reads(controlcard,"FOURIER",fouriername_t,fouriername)
183 open (ifourier,file=fouriername_t,status='old')
185 call reads(controlcard,"ELEPAR",elename_t,elename)
186 open (ielep,file=elename_t,status='old')
188 call reads(controlcard,"SIDEPAR",sidename_t,sidename)
189 open (isidep,file=sidename_t,status='old')
191 call reads(controlcard,"SCPPAR",scpname_t,scpname)
192 open (iscpp,file=scpname_t,status='old')
194 write (iout,*) "Parameter set:",iparm
195 write (iout,*) "Energy-term weights:"
197 write (iout,'(a16,f10.5)') wname(i),ww(i)
199 write (iout,*) "Sidechain potential file : ",
200 & sidename_t(:ilen(sidename_t))
202 write (iout,*) "SCp potential file : ",
203 & scpname_t(:ilen(scpname_t))
205 write (iout,*) "Electrostatic potential file : ",
206 & elename_t(:ilen(elename_t))
207 write (iout,*) "Cumulant coefficient file : ",
208 & fouriername_t(:ilen(fouriername_t))
209 write (iout,*) "Torsional parameter file : ",
210 & torname_t(:ilen(torname_t))
211 write (iout,*) "Double torsional parameter file : ",
212 & tordname_t(:ilen(tordname_t))
213 write (iout,*) "Backbone-rotamer parameter file : ",
214 & sccorname(:ilen(sccorname))
215 write (iout,*) "Bond & inertia constant file : ",
216 & bondname_t(:ilen(bondname_t))
217 write (iout,*) "Bending parameter file : ",
218 & thetname_t(:ilen(thetname_t))
219 write (iout,*) "Rotamer parameter file : ",
220 & rotname_t(:ilen(rotname_t))
223 c Read the virtual-bond parameters, masses, and moments of inertia
224 c and Stokes' radii of the peptide group and side chains
227 read (ibond,*,end=121,err=121) vbldp0,vbldpdum,akp
230 read (ibond,*,end=121,err=121) vbldsc0(1,i),aksc(1,i)
231 dsc(i) = vbldsc0(1,i)
235 dsc_inv(i)=1.0D0/dsc(i)
239 read (ibond,*,end=121,err=121) ijunk,vbldp0,vbldpdum,akp,rjunk
241 read (ibond,*,end=121,err=121) nbondterm(i),(vbldsc0(j,i),
242 & aksc(j,i),abond0(j,i),j=1,nbondterm(i))
243 dsc(i) = vbldsc0(1,i)
247 dsc_inv(i)=1.0D0/dsc(i)
252 write(iout,'(/a/)')"Force constants virtual bonds:"
253 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
255 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
257 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
258 & vbldsc0(1,i),aksc(1,i),abond0(1,i)
260 write (iout,'(13x,3f10.5)')
261 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
265 read(iliptranpar,*,end=1161,err=1161) pepliptran
267 read(iliptranpar,*,end=1161,err=1161) liptranene(i)
272 C Read the parameters of the probability distribution/energy expression
273 C of the virtual-bond valence angles theta
276 read (ithep,*,end=111,err=111) a0thet(i),(athet(j,i,1,1),j=1,2),
277 & (bthet(j,i,1,1),j=1,2)
278 read (ithep,*,end=111,err=111) (polthet(j,i),j=0,3)
279 read (ithep,*,end=111,err=111) (gthet(j,i),j=1,3)
280 read (ithep,*,end=111,err=111) theta0(i),sig0(i),sigc0(i)
284 athet(1,i,1,-1)=athet(1,i,1,1)
285 athet(2,i,1,-1)=athet(2,i,1,1)
286 bthet(1,i,1,-1)=-bthet(1,i,1,1)
287 bthet(2,i,1,-1)=-bthet(2,i,1,1)
288 athet(1,i,-1,1)=-athet(1,i,1,1)
289 athet(2,i,-1,1)=-athet(2,i,1,1)
290 bthet(1,i,-1,1)=bthet(1,i,1,1)
291 bthet(2,i,-1,1)=bthet(2,i,1,1)
295 athet(1,i,-1,-1)=athet(1,-i,1,1)
296 athet(2,i,-1,-1)=-athet(2,-i,1,1)
297 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
298 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
299 athet(1,i,-1,1)=athet(1,-i,1,1)
300 athet(2,i,-1,1)=-athet(2,-i,1,1)
301 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
302 bthet(2,i,-1,1)=bthet(2,-i,1,1)
303 athet(1,i,1,-1)=-athet(1,-i,1,1)
304 athet(2,i,1,-1)=athet(2,-i,1,1)
305 bthet(1,i,1,-1)=bthet(1,-i,1,1)
306 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
311 polthet(j,i)=polthet(j,-i)
314 gthet(j,i)=gthet(j,-i)
320 c & 'Parameters of the virtual-bond valence angles:'
321 c write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
322 c & ' ATHETA0 ',' A1 ',' A2 ',
325 c write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
326 c & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
328 c write (iout,'(/a/9x,5a/79(1h-))')
329 c & 'Parameters of the expression for sigma(theta_c):',
330 c & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
331 c & ' ALPH3 ',' SIGMA0C '
333 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
334 c & (polthet(j,i),j=0,3),sigc0(i)
336 c write (iout,'(/a/9x,5a/79(1h-))')
337 c & 'Parameters of the second gaussian:',
338 c & ' THETA0 ',' SIGMA0 ',' G1 ',
341 c write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
342 c & sig0(i),(gthet(j,i),j=1,3)
345 & 'Parameters of the virtual-bond valence angles:'
346 write (iout,'(/a/9x,5a/79(1h-))')
347 & 'Coefficients of expansion',
348 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
349 & ' b1*10^1 ',' b2*10^1 '
351 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
352 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
353 & (10*bthet(j,i,1,1),j=1,2)
355 write (iout,'(/a/9x,5a/79(1h-))')
356 & 'Parameters of the expression for sigma(theta_c):',
357 & ' alpha0 ',' alph1 ',' alph2 ',
358 & ' alhp3 ',' sigma0c '
360 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
361 & (polthet(j,i),j=0,3),sigc0(i)
363 write (iout,'(/a/9x,5a/79(1h-))')
364 & 'Parameters of the second gaussian:',
365 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
368 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
369 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
373 IF (tor_mode.eq.0) THEN
375 C Read the parameters of Utheta determined from ab initio surfaces
376 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
378 read (ithep,*,end=111,err=111) nthetyp,ntheterm,ntheterm2,
379 & ntheterm3,nsingle,ndouble
380 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
381 read (ithep,*,end=111,err=111) (ithetyp(i),i=1,ntyp1)
383 ithetyp(i)=-ithetyp(-i)
386 do i=-maxthetyp,maxthetyp
387 do j=-maxthetyp,maxthetyp
388 do k=-maxthetyp,maxthetyp
389 aa0thet(i,j,k,iblock)=0.0d0
391 aathet(l,i,j,k,iblock)=0.0d0
395 bbthet(m,l,i,j,k,iblock)=0.0d0
396 ccthet(m,l,i,j,k,iblock)=0.0d0
397 ddthet(m,l,i,j,k,iblock)=0.0d0
398 eethet(m,l,i,j,k,iblock)=0.0d0
404 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
405 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
413 C write (iout,*) "KURWA1"
416 do j=-nthetyp,nthetyp
417 do k=-nthetyp,nthetyp
418 read (ithep,'(6a)',end=111,err=111) res1
419 write(iout,*) res1,i,j,k
420 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
421 read (ithep,*,end=111,err=111)(aathet(l,i,j,k,iblock),
423 read (ithep,*,end=111,err=111)
424 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
425 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
426 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
427 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle)
429 read (ithep,*,end=111,err=111)
430 & (((ffthet(llll,lll,ll,i,j,k,iblock),
431 & ffthet(lll,llll,ll,i,j,k,iblock),
432 & ggthet(llll,lll,ll,i,j,k,iblock)
433 & ,ggthet(lll,llll,ll,i,j,k,iblock),
434 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
438 C write(iout,*) "KURWA1.1"
440 C For dummy ends assign glycine-type coefficients of theta-only terms; the
441 C coefficients of theta-and-gamma-dependent terms are zero.
446 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
447 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
449 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
450 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
453 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
455 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
458 C write(iout,*) "KURWA1.5"
459 C Substitution for D aminoacids from symmetry.
462 do j=-nthetyp,nthetyp
463 do k=-nthetyp,nthetyp
464 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
466 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
470 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
471 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
472 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
473 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
479 ffthet(llll,lll,ll,i,j,k,iblock)=
480 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
481 ffthet(lll,llll,ll,i,j,k,iblock)=
482 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
483 ggthet(llll,lll,ll,i,j,k,iblock)=
484 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
485 ggthet(lll,llll,ll,i,j,k,iblock)=
486 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
496 C Control printout of the coefficients of virtual-bond-angle potentials
499 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
504 write (iout,'(//4a)')
505 & 'Type ',onelett(i),onelett(j),onelett(k)
506 write (iout,'(//a,10x,a)') " l","a[l]"
507 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
508 write (iout,'(i2,1pe15.5)')
509 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
511 write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
512 & "b",l,"c",l,"d",l,"e",l
514 write (iout,'(i2,4(1pe15.5))') m,
515 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
516 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
520 write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
521 & "f+",l,"f-",l,"g+",l,"g-",l
524 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
525 & ffthet(n,m,l,i,j,k,iblock),
526 & ffthet(m,n,l,i,j,k,iblock),
527 & ggthet(n,m,l,i,j,k,iblock),
528 & ggthet(m,n,l,i,j,k,iblock)
541 C here will be the apropriate recalibrating for D-aminoacid
542 read (ithep,*,end=111,err=111) nthetyp
543 do i=-nthetyp+1,nthetyp-1
544 read (ithep,*,end=111,err=111) nbend_kcc_Tb(i)
545 do j=0,nbend_kcc_Tb(i)
546 read (ithep,*,end=111,err=111) ijunk,v1bend_chyb(j,i)
551 & "Parameters of the valence-only potentials"
552 do i=-nthetyp+1,nthetyp-1
553 write (iout,'(2a)') "Type ",toronelet(i)
554 do k=0,nbend_kcc_Tb(i)
555 write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
564 C write(iout,*) 'KURWA2'
567 C Read the parameters of the probability distribution/energy expression
568 C of the side chains.
571 cc write (iout,*) "tu dochodze",i
572 read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
576 dsc_inv(i)=1.0D0/dsc(i)
587 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
588 & ((blower(k,l,1),l=1,k),k=1,3)
589 censc(1,1,-i)=censc(1,1,i)
590 censc(2,1,-i)=censc(2,1,i)
591 censc(3,1,-i)=-censc(3,1,i)
593 read (irotam,*,end=112,err=112) bsc(j,i)
594 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
595 & ((blower(k,l,j),l=1,k),k=1,3)
596 censc(1,j,-i)=censc(1,j,i)
597 censc(2,j,-i)=censc(2,j,i)
598 censc(3,j,-i)=-censc(3,j,i)
599 C BSC is amplitude of Gaussian
606 akl=akl+blower(k,m,j)*blower(l,m,j)
610 if (((k.eq.3).and.(l.ne.3))
611 & .or.((l.eq.3).and.(k.ne.3))) then
612 gaussc(k,l,j,-i)=-akl
613 gaussc(l,k,j,-i)=-akl
625 write (iout,'(/a)') 'Parameters of side-chain local geometry'
629 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
630 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
631 c write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
632 c write (iout,'(a,f10.4,4(16x,f10.4))')
633 c & 'Center ',(bsc(j,i),j=1,nlobi)
634 c write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
635 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
636 & 'log h',(bsc(j,i),j=1,nlobi)
637 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
638 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
645 c blower(k,l,j)=gaussc(ind,j,i)
650 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
651 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
658 C Read scrot parameters for potentials determined from all-atom AM1 calculations
659 C added by Urszula Kozlowska 07/11/2007
662 read (irotam,*,end=112,err=112)
664 read (irotam,*,end=112,err=112)
667 read(irotam,*,end=112,err=112) sc_parmin(j,i)
674 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
675 C interaction energy of the Gly, Ala, and Pro prototypes.
677 read (ifourier,*,end=115,err=115) nloctyp
679 read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
680 read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
681 itype2loc(ntyp1)=nloctyp
682 iloctyp(nloctyp)=ntyp1
684 itype2loc(-i)=-itype2loc(i)
693 iloctyp(-i)=-iloctyp(i)
695 c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
696 c write (iout,*) "nloctyp",nloctyp,
697 c & " iloctyp",(iloctyp(i),i=0,nloctyp)
700 c write (iout,*) "NEWCORR",i
701 read (ifourier,*,end=115,err=115)
704 read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
707 c write (iout,*) "NEWCORR BNEW1"
708 c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
711 read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
714 c write (iout,*) "NEWCORR BNEW2"
715 c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
717 read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
718 read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
720 c write (iout,*) "NEWCORR CCNEW"
721 c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
723 read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
724 read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
726 c write (iout,*) "NEWCORR DDNEW"
727 c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
731 read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
735 c write (iout,*) "NEWCORR EENEW1"
736 c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
738 read (ifourier,*,end=115,err=115) e0new(ii,i)
740 c write (iout,*) (e0new(ii,i),ii=1,3)
742 c write (iout,*) "NEWCORR EENEW"
745 ccnew(ii,1,i)=ccnew(ii,1,i)/2
746 ccnew(ii,2,i)=ccnew(ii,2,i)/2
747 ddnew(ii,1,i)=ddnew(ii,1,i)/2
748 ddnew(ii,2,i)=ddnew(ii,2,i)/2
753 bnew1(ii,1,-i)= bnew1(ii,1,i)
754 bnew1(ii,2,-i)=-bnew1(ii,2,i)
755 bnew2(ii,1,-i)= bnew2(ii,1,i)
756 bnew2(ii,2,-i)=-bnew2(ii,2,i)
759 c ccnew(ii,1,i)=ccnew(ii,1,i)/2
760 c ccnew(ii,2,i)=ccnew(ii,2,i)/2
761 c ddnew(ii,1,i)=ddnew(ii,1,i)/2
762 c ddnew(ii,2,i)=ddnew(ii,2,i)/2
763 ccnew(ii,1,-i)=ccnew(ii,1,i)
764 ccnew(ii,2,-i)=-ccnew(ii,2,i)
765 ddnew(ii,1,-i)=ddnew(ii,1,i)
766 ddnew(ii,2,-i)=-ddnew(ii,2,i)
768 e0new(1,-i)= e0new(1,i)
769 e0new(2,-i)=-e0new(2,i)
770 e0new(3,-i)=-e0new(3,i)
772 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
773 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
774 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
775 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
779 write (iout,'(a)') "Coefficients of the multibody terms"
780 do i=-nloctyp+1,nloctyp-1
781 write (iout,*) "Type: ",onelet(iloctyp(i))
782 write (iout,*) "Coefficients of the expansion of B1"
784 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
786 write (iout,*) "Coefficients of the expansion of B2"
788 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
790 write (iout,*) "Coefficients of the expansion of C"
791 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
792 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
793 write (iout,*) "Coefficients of the expansion of D"
794 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
795 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
796 write (iout,*) "Coefficients of the expansion of E"
797 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
800 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
807 & write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
809 read (ifourier,*,end=115,err=115)
810 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
812 write (iout,*) 'Type ',onelet(iloctyp(i))
813 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
827 c B1tilde(1,i) = b(3)
828 c B1tilde(2,i) =-b(5)
829 c B1tilde(1,-i) =-b(3)
830 c B1tilde(2,-i) =b(5)
837 cc B1tilde(1,i) = b(3,i)
838 cc B1tilde(2,i) =-b(5,i)
839 C B1tilde(1,-i) =-b(3,i)
840 C B1tilde(2,-i) =b(5,i)
841 cc b1tilde(1,i)=0.0d0
842 cc b1tilde(2,i)=0.0d0
854 CCold(1,1,-i)= b(7,i)
855 CCold(2,2,-i)=-b(7,i)
856 CCold(2,1,-i)=-b(9,i)
857 CCold(1,2,-i)=-b(9,i)
862 c Ctilde(1,1,i)= CCold(1,1,i)
863 c Ctilde(1,2,i)= CCold(1,2,i)
864 c Ctilde(2,1,i)=-CCold(2,1,i)
865 c Ctilde(2,2,i)=-CCold(2,2,i)
867 c Ctilde(1,1,i)=0.0d0
868 c Ctilde(1,2,i)=0.0d0
869 c Ctilde(2,1,i)=0.0d0
870 c Ctilde(2,2,i)=0.0d0
875 DDold(1,1,-i)= b(6,i)
876 DDold(2,2,-i)=-b(6,i)
877 DDold(2,1,-i)=-b(8,i)
878 DDold(1,2,-i)=-b(8,i)
883 c Dtilde(1,1,i)= DD(1,1,i)
884 c Dtilde(1,2,i)= DD(1,2,i)
885 c Dtilde(2,1,i)=-DD(2,1,i)
886 c Dtilde(2,2,i)=-DD(2,2,i)
888 c Dtilde(1,1,i)=0.0d0
889 c Dtilde(1,2,i)=0.0d0
890 c Dtilde(2,1,i)=0.0d0
891 c Dtilde(2,2,i)=0.0d0
892 EEold(1,1,i)= b(10,i)+b(11,i)
893 EEold(2,2,i)=-b(10,i)+b(11,i)
894 EEold(2,1,i)= b(12,i)-b(13,i)
895 EEold(1,2,i)= b(12,i)+b(13,i)
896 EEold(1,1,-i)= b(10,i)+b(11,i)
897 EEold(2,2,-i)=-b(10,i)+b(11,i)
898 EEold(2,1,-i)=-b(12,i)+b(13,i)
899 EEold(1,2,-i)=-b(12,i)-b(13,i)
900 write(iout,*) "TU DOCHODZE"
906 c ee(2,1,i)=ee(1,2,i)
911 &"Coefficients of the cumulants (independent of valence angles)"
912 do i=-nloctyp+1,nloctyp-1
913 write (iout,*) 'Type ',onelet(iloctyp(i))
915 write(iout,'(2f10.5)') B(3,i),B(5,i)
917 write(iout,'(2f10.5)') B(2,i),B(4,i)
920 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
924 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
928 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
933 C write (iout,*) 'KURWAKURWA'
936 C Read torsional parameters in old format
938 read (itorp,*,end=113,err=113) ntortyp,nterm_old
939 write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
940 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
943 read (itorp,'(a)',end=113,err=113)
945 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
951 write (iout,'(/a/)') 'Torsional constants:'
954 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
955 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
961 C Read torsional parameters
963 IF (TOR_MODE.eq.0) THEN
965 read (itorp,*,end=113,err=113) ntortyp
966 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
967 write (iout,*) 'ntortyp',ntortyp
969 itype2loc(i)=itortyp(i)
972 itype2loc(-i)=-itype2loc(i)
976 itortyp(i)=-itortyp(-i)
978 c write (iout,*) 'ntortyp',ntortyp
980 do j=-ntortyp+1,ntortyp-1
981 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
983 nterm(-i,-j,iblock)=nterm(i,j,iblock)
984 nlor(-i,-j,iblock)=nlor(i,j,iblock)
987 do k=1,nterm(i,j,iblock)
988 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
990 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
991 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
992 v0ij=v0ij+si*v1(k,i,j,iblock)
995 do k=1,nlor(i,j,iblock)
996 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
997 & vlor2(k,i,j),vlor3(k,i,j)
998 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1001 v0(-i,-j,iblock)=v0ij
1007 write (iout,'(/a/)') 'Torsional constants:'
1011 write (iout,*) 'ityp',i,' jtyp',j
1012 write (iout,*) 'Fourier constants'
1013 do k=1,nterm(i,j,iblock)
1014 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
1017 write (iout,*) 'Lorenz constants'
1018 do k=1,nlor(i,j,iblock)
1019 write (iout,'(3(1pe15.5))')
1020 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1027 C 6/23/01 Read parameters for double torsionals
1031 do j=-ntortyp+1,ntortyp-1
1032 do k=-ntortyp+1,ntortyp-1
1033 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1034 c write (iout,*) "OK onelett",
1037 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
1038 & .or. t3.ne.toronelet(k)) then
1039 write (iout,*) "Error in double torsional parameter file",
1042 call MPI_Finalize(Ierror)
1044 stop "Error in double torsional parameter file"
1046 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
1047 & ntermd_2(i,j,k,iblock)
1048 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1049 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1050 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
1051 & ntermd_1(i,j,k,iblock))
1052 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
1053 & ntermd_1(i,j,k,iblock))
1054 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
1055 & ntermd_1(i,j,k,iblock))
1056 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
1057 & ntermd_1(i,j,k,iblock))
1058 C Martix of D parameters for one dimesional foureir series
1059 do l=1,ntermd_1(i,j,k,iblock)
1060 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1061 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1062 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1063 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1064 c write(iout,*) "whcodze" ,
1065 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1067 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
1068 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
1069 & v2s(m,l,i,j,k,iblock),
1070 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1071 C Martix of D parameters for two dimesional fourier series
1072 do l=1,ntermd_2(i,j,k,iblock)
1074 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1075 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1076 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1077 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1086 write (iout,*) 'Constants for double torsionals'
1089 do j=-ntortyp+1,ntortyp-1
1090 do k=-ntortyp+1,ntortyp-1
1091 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
1092 & ' nsingle',ntermd_1(i,j,k,iblock),
1093 & ' ndouble',ntermd_2(i,j,k,iblock)
1095 write (iout,*) 'Single angles:'
1096 do l=1,ntermd_1(i,j,k,iblock)
1097 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
1098 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
1099 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
1100 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1103 write (iout,*) 'Pairs of angles:'
1104 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1105 do l=1,ntermd_2(i,j,k,iblock)
1106 write (iout,'(i5,20f10.5)')
1107 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1110 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1111 do l=1,ntermd_2(i,j,k,iblock)
1112 write (iout,'(i5,20f10.5)')
1113 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
1114 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1123 ELSE IF (TOR_MODE.eq.1) THEN
1125 C read valence-torsional parameters
1126 read (itorp,*,end=113,err=113) ntortyp
1128 write (iout,*) "Valence-torsional parameters read in ntortyp",
1130 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1131 write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
1133 itortyp(i)=-itortyp(-i)
1135 do i=-ntortyp+1,ntortyp-1
1136 do j=-ntortyp+1,ntortyp-1
1137 C first we read the cos and sin gamma parameters
1138 read (itorp,'(13x,a)',end=113,err=113) string
1139 write (iout,*) i,j,string
1140 read (itorp,*,end=113,err=113)
1141 & nterm_kcc(j,i),nterm_kcc_Tb(j,i)
1142 C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
1143 do k=1,nterm_kcc(j,i)
1144 do l=1,nterm_kcc_Tb(j,i)
1145 do ll=1,nterm_kcc_Tb(j,i)
1146 read (itorp,*,end=113,err=113) ii,jj,kk,
1147 & v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1155 c AL 4/8/16: Calculate coefficients from one-body parameters
1158 itortyp(i)=itype2loc(i)
1162 &"Val-tor parameters calculated from cumulant coefficients ntortyp"
1164 do i=-ntortyp+1,ntortyp-1
1165 do j=-ntortyp+1,ntortyp-1
1168 do k=1,nterm_kcc_Tb(j,i)
1169 do l=1,nterm_kcc_Tb(j,i)
1170 v1_kcc(k,l,1,i,j)=bnew1(k,1,i)*bnew2(l,1,j)
1171 & +bnew1(k,2,i)*bnew2(l,2,j)
1172 v2_kcc(k,l,1,i,j)=bnew1(k,1,i)*bnew2(l,2,j)
1173 & +bnew1(k,2,i)*bnew2(l,1,j)
1176 do k=1,nterm_kcc_Tb(j,i)
1177 do l=1,nterm_kcc_Tb(j,i)
1178 v1_kcc(k,l,2,i,j)=-0.25*(ccnew(k,1,i)*ddnew(l,1,j)
1179 & -ccnew(k,2,i)*ddnew(l,2,j))
1180 v2_kcc(k,l,2,i,j)=-0.25*(ccnew(k,2,i)*ddnew(l,1,j)
1181 & +ccnew(k,1,i)*ddnew(l,2,j))
1184 cf(theta,gamma)=-(b21(theta)*b11(theta)+b12(theta)*b22(theta))*cos(gamma)-(b22(theta)*b11(theta)+b21(theta)*b12(theta))*sin(gamma)+(c11(theta)*d11(theta)-c12(theta)*d12(theta))*cos(2*gamma)+(c12(theta)*d11(theta)+c11(theta)*d12(theta))*sin(2*gamma)
1190 if (tor_mode.gt.0 .and. lprint) then
1191 c Print valence-torsional parameters
1193 & "Parameters of the valence-torsional potentials"
1194 do i=-ntortyp+1,ntortyp-1
1195 do j=-ntortyp+1,ntortyp-1
1196 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
1197 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
1198 do k=1,nterm_kcc(j,i)
1199 do l=1,nterm_kcc_Tb(j,i)
1200 do ll=1,nterm_kcc_Tb(j,i)
1201 write (iout,'(3i5,2f15.4)')
1202 & k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1211 C Read of Side-chain backbone correlation parameters
1212 C Modified 11 May 2012 by Adasko
1215 read (isccor,*,end=119,err=119) nsccortyp
1216 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1218 isccortyp(i)=-isccortyp(-i)
1220 iscprol=isccortyp(20)
1221 c write (iout,*) 'ntortyp',ntortyp
1223 cc maxinter is maximum interaction sites
1227 read (isccor,*,end=119,err=119)
1228 &nterm_sccor(i,j),nlor_sccor(i,j)
1229 c write (iout,*) nterm_sccor(i,j)
1235 nterm_sccor(-i,j)=nterm_sccor(i,j)
1236 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1237 nterm_sccor(i,-j)=nterm_sccor(i,j)
1238 c write (iout,*) nterm_sccor(i,j),nterm_sccor(-i,j),
1239 c & nterm_sccor(-i,-j),nterm_sccor(i,-j)
1240 do k=1,nterm_sccor(i,j)
1241 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
1243 if (j.eq.iscprol) then
1244 if (i.eq.isccortyp(10)) then
1245 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1246 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1248 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
1249 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1250 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
1251 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1252 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1253 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1254 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1255 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1258 if (i.eq.isccortyp(10)) then
1259 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1260 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1262 if (j.eq.isccortyp(10)) then
1263 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1264 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1266 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1267 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1268 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1269 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1270 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1271 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1275 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1276 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1277 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1278 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1281 do k=1,nlor_sccor(i,j)
1282 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
1283 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1284 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
1285 &(1+vlor3sccor(k,i,j)**2)
1287 v0sccor(l,i,j)=v0ijsccor
1288 v0sccor(l,-i,j)=v0ijsccor1
1289 v0sccor(l,i,-j)=v0ijsccor2
1290 v0sccor(l,-i,-j)=v0ijsccor3
1296 write (iout,'(/a/)') 'Torsional constants of SCCORR:'
1300 write (iout,*) 'ityp',i,' jtyp',j
1301 write (iout,*) 'Fourier constants'
1302 do k=1,nterm_sccor(i,j)
1303 write (iout,'(2(1pe15.5))')
1304 & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1306 write (iout,*) 'Lorenz constants'
1307 do k=1,nlor_sccor(i,j)
1308 write (iout,'(3(1pe15.5))')
1309 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1316 C Read electrostatic-interaction parameters
1319 write (iout,'(/a)') 'Electrostatic interaction constants:'
1320 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1321 & 'IT','JT','APP','BPP','AEL6','AEL3'
1323 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1324 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1325 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1326 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1331 app (i,j)=epp(i,j)*rri*rri
1332 bpp (i,j)=-2.0D0*epp(i,j)*rri
1333 ael6(i,j)=elpp6(i,j)*4.2D0**6
1334 ael3(i,j)=elpp3(i,j)*4.2D0**3
1336 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1337 & ael6(i,j),ael3(i,j)
1341 C Read side-chain interaction parameters.
1343 read (isidep,*) ipot,expon
1344 if (ipot.lt.1 .or. ipot.gt.6) then
1345 write (iout,'(2a)') 'Error while reading SC interaction',
1346 & 'potential file - unknown potential type.'
1350 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1351 & ', exponents are ',expon,2*expon
1352 goto (10,20,30,30,40,50,60) ipot
1353 C----------------------- LJ potential ---------------------------------
1354 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1356 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1357 write (iout,'(a/)') 'The epsilon array:'
1358 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1359 write (iout,'(/a)') 'One-body parameters:'
1360 write (iout,'(a,4x,a)') 'residue','sigma'
1361 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1364 C----------------------- LJK potential --------------------------------
1365 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1366 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1368 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1369 write (iout,'(a/)') 'The epsilon array:'
1370 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1371 write (iout,'(/a)') 'One-body parameters:'
1372 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1373 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1377 C---------------------- GB or BP potential -----------------------------
1378 30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1379 & (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
1381 C For the GB potential convert sigma'**2 into chi'
1384 chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
1388 write (iout,'(/a/)') 'Parameters of the BP potential:'
1389 write (iout,'(a/)') 'The epsilon array:'
1390 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1391 write (iout,'(/a)') 'One-body parameters:'
1392 write (iout,'(a,4x,5a)') 'residue',' sigma ','s||/s_|_^2',
1393 & ' chip0 ',' chip ',' alph '
1394 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),sigii(i),
1395 & chip0(i),chip(i),alp(i),i=1,ntyp)
1398 C--------------------- GBV potential -----------------------------------
1399 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1400 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1401 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1403 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1404 write (iout,'(a/)') 'The epsilon array:'
1405 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1406 write (iout,'(/a)') 'One-body parameters:'
1407 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1408 & 's||/s_|_^2',' chip ',' alph '
1409 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1410 & sigii(i),chip(i),alp(i),i=1,ntyp)
1413 C--------------------- Momo potential -----------------------------------
1415 read (isidep,*) (icharge(i),i=1,ntyp)
1416 c write (2,*) "icharge",(icharge(i),i=1,ntyp)
1419 c! write (*,*) "Im in ", i, " ", j
1421 & eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i),
1422 & (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j),
1423 & chis(i,j),chis(j,i),
1424 & nstate(i,j),(wstate(k,i,j),k=1,4),
1429 & dtail(1,i,j),dtail(2,i,j),
1430 & epshead(i,j),sig0head(i,j),
1431 & rborn(i,j),rborn(j,i),
1432 & (wqdip(k,i,j),k=1,2),wquad(i,j),
1433 & alphapol(i,j),alphapol(j,i),
1434 & (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j)
1437 c! write (*,*) "nstate gly-gly", nstate(10,10)
1438 c! THIS LOOP FILLS PARAMETERS FOR PAIRS OF AA's NOT EXPLICITLY
1439 c! DEFINED IN SCPARM.MOMO. IT DOES SO BY TAKING THEM FROM SYMMETRIC
1440 c! PAIR, FOR EG. IF ARG-HIS IS BLANK, IT WILL TAKE PARAMETERS
1443 c! DISABLE IT AT >>YOUR OWN PERIL<<
1448 sigma(i,j) = sigma(j,i)
1449 nstate(i,j) = nstate(j,i)
1450 sigmap1(i,j) = sigmap1(j,i)
1451 sigmap2(i,j) = sigmap2(j,i)
1452 sigiso1(i,j) = sigiso1(j,i)
1453 sigiso2(i,j) = sigiso2(j,i)
1456 alphasur(k,i,j) = alphasur(k,j,i)
1457 wstate(k,i,j) = wstate(k,j,i)
1458 alphiso(k,i,j) = alphiso(k,j,i)
1461 dhead(2,1,i,j) = dhead(1,1,j,i)
1462 dhead(2,2,i,j) = dhead(1,2,j,i)
1463 dhead(1,1,i,j) = dhead(2,1,j,i)
1464 dhead(1,2,i,j) = dhead(2,2,j,i)
1465 dtail(1,i,j) = dtail(1,j,i)
1466 dtail(2,i,j) = dtail(2,j,i)
1469 c! write (*,*) "dhead(k,l,j,i) = ", dhead(k,l,j,i)
1470 c! write (*,*) "dhead(k,l,i,j) = ", dhead(k,l,i,j)
1471 c! dhead(l,k,i,j) = dhead(k,l,j,i)
1475 epshead(i,j) = epshead(j,i)
1476 sig0head(i,j) = sig0head(j,i)
1479 wqdip(k,i,j) = wqdip(k,j,i)
1482 wquad(i,j) = wquad(j,i)
1483 epsintab(i,j) = epsintab(j,i)
1488 if (.not.lprint) goto 70
1490 & "Parameters of the new physics-based SC-SC interaction potential"
1491 write (iout,'(/7a)') 'Residues',' epsGB',' rGB',
1492 & ' chi1GB',' chi2GB',' chip1GB',' chip2GB'
1495 write (iout,'(2(a3,1x),1pe10.3,5(0pf10.3))')
1496 & restyp(i),restyp(j),eps(i,j),sigma(i,j),chi(i,j),chi(j,i),
1497 & chipp(i,j),chipp(j,i)
1500 write (iout,'(/9a)') 'Residues',' alphasur1',' alphasur2',
1501 & ' alphasur3',' alphasur4',' sigmap1',' sigmap2',
1505 write (iout,'(2(a3,1x),8(0pf10.3))')
1506 & restyp(i),restyp(j),(alphasur(k,i,j),k=1,4),
1507 & sigmap1(i,j),sigmap2(j,i),chis(i,j),chis(j,i)
1510 write (iout,'(/14a)') 'Residues',' nst ',' wst1',
1511 & ' wst2',' wst3',' wst4',' dh11',' dh21',
1512 & ' dh12',' dh22',' dt1',' dt2',' epsh1',
1516 write (iout,'(2(a3,1x),i3,4f8.3,6f7.2,f9.5,f7.2)')
1517 & restyp(i),restyp(j),nstate(i,j),(wstate(k,i,j),k=1,4),
1518 & ((dhead(l,k,i,j),l=1,2),k=1,2),dtail(1,i,j),dtail(2,i,j),
1519 & epshead(i,j),sig0head(i,j)
1522 write (iout,'(/12a)') 'Residues',' ch1',' ch2',
1523 & ' rborn1',' rborn2',' wqdip1',' wqdip2',
1527 write (iout,'(2(a3,1x),2i4,5f10.3)')
1528 & restyp(i),restyp(j),icharge(i),icharge(j),
1529 & rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j)
1532 write (iout,'(/12a)') 'Residues',
1534 & ' alphpol2',' alphiso1',' alpiso2',
1535 & ' alpiso3',' alpiso4',' sigiso1',' sigiso2',
1539 write (iout,'(2(a3,1x),11f10.3)')
1540 & restyp(i),restyp(j),alphapol(i,j),alphapol(j,i),
1541 & (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(j,i),
1546 C For the GB potential convert sigma'**2 into chi'
1548 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1551 write (iout,'(/a/)') 'Parameters of the BP potential:'
1552 write (iout,'(a/)') 'The epsilon array:'
1553 CALL printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1554 write (iout,'(/a)') 'One-body parameters:'
1555 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1557 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1558 & chip(i),alp(i),i=1,ntyp)
1562 C-----------------------------------------------------------------------
1563 C Calculate the "working" parameters of SC interactions.
1568 alphasur(k,j,i)=alphasur(k,i,j)
1569 alphiso(k,j,i)=alphiso(k,i,j)
1570 wstate(k,j,i)=wstate(k,i,j)
1573 wqdip(k,j,i)=wqdip(k,i,j)
1577 dhead(l,k,j,i)=dhead(l,k,i,j)
1585 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1586 sigma(j,i)=sigma(i,j)
1587 rs0(i,j)=dwa16*sigma(i,j)
1594 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1595 & 'Working parameters of the SC interactions:',
1596 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1601 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4 .or. ipot.eq.6 )
1611 sigeps=dsign(1.0D0,epsij)
1613 aa(i,j)=epsij*rrij*rrij
1614 bb(i,j)=-sigeps*epsij*rrij
1617 IF ((ipot.gt.2).AND.(ipot.LT.6)) THEN
1618 sigt1sq=sigma0(i)**2
1619 sigt2sq=sigma0(j)**2
1622 ratsig1=sigt2sq/sigt1sq
1623 ratsig2=1.0D0/ratsig1
1624 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1625 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1626 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1630 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1631 sigmaii(i,j)=rsum_max
1632 sigmaii(j,i)=rsum_max
1634 c sigmaii(i,j)=r0(i,j)
1635 c sigmaii(j,i)=r0(i,j)
1637 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1638 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1639 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1640 augm(i,j)=epsij*r_augm**(2*expon)
1641 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1649 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1650 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1651 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1653 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3),2i3,10f8.4,
1655 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1656 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i),
1657 & icharge(i),icharge(j),chipp(i,j),chipp(j,i),
1658 & (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(j,i),
1659 & chis(i,j),chis(j,i),
1660 & nstate(i,j),(wstate(k,i,j),k=1,4),
1661 & ((dhead(l,k,i,j),l=1,2),k=1,2),dtail(1,i,j),dtail(2,i,j),
1662 & epshead(i,j),sig0head(i,j),
1663 & rborn(i,j),(wqdip(k,i,j),k=1,2),wquad(i,j),
1664 & alphapol(i,j),alphapol(j,i),
1665 & (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j)
1672 C Define the SC-p interaction constants
1676 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1678 c aad(i,1)=0.3D0*4.0D0**12
1679 C Following line for constants currently implemented
1680 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1681 aad(i,1)=1.5D0*4.0D0**12
1682 c aad(i,1)=0.17D0*5.6D0**12
1684 C "Soft" SC-p repulsion
1686 C Following line for constants currently implemented
1687 c aad(i,1)=0.3D0*4.0D0**6
1688 C "Hard" SC-p repulsion
1689 bad(i,1)=3.0D0*4.0D0**6
1690 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1699 C 8/9/01 Read the SC-p interaction constants from file
1702 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1705 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1706 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1707 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1708 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1712 write (iout,*) "Parameters of SC-p interactions:"
1714 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1715 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1720 C Define the constants of the disulfide bridge
1724 c Old arbitrary potential - commented out.
1729 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1730 c energy surface of diethyl disulfide.
1731 c A. Liwo and U. Kozlowska, 11/24/03
1740 write (iout,*) dyn_ss,'dyndyn'
1742 ss_depth=ebr/wsc-0.25*eps(1,1)
1743 C write(iout,*) akcm,whpb,wsc,'KURWA'
1744 Ht=Ht/wsc-0.25*eps(1,1)
1753 ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
1757 write (iout,'(/a)') "Disulfide bridge parameters:"
1758 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1759 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1760 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1761 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1764 if (shield_mode.gt.0) then
1766 C VSolvSphere the volume of solving sphere
1768 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1769 C there will be no distinction between proline peptide group and normal peptide
1770 C group in case of shielding parameters
1771 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1772 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1773 write (iout,*) VSolvSphere,VSolvSphere_div
1774 C long axis of side chain
1776 long_r_sidechain(i)=vbldsc0(1,i)
1777 short_r_sidechain(i)=sigma0(i)
1782 111 write (iout,*) "Error reading bending energy parameters."
1784 112 write (iout,*) "Error reading rotamer energy parameters."
1786 113 write (iout,*) "Error reading torsional energy parameters."
1788 114 write (iout,*) "Error reading double torsional energy parameters."
1791 & "Error reading cumulant (multibody energy) parameters."
1793 116 write (iout,*) "Error reading electrostatic energy parameters."
1795 1161 write (iout,*) "Error reading lipid energy parameters."
1797 117 write (iout,*) "Error reading side chain interaction parameters."
1799 118 write (iout,*) "Error reading SCp interaction parameters."
1801 119 write (iout,*) "Error reading SCCOR parameters"
1803 121 write (iout,*) "Error reading bond parameters"
1806 call MPI_Finalize(Ierror)