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 C Important! Energy-term weights ARE NOT read here; they are read from the
7 C main input file instead, because NO defaults have yet been set for these
16 include 'COMMON.IOUNITS'
17 include 'COMMON.CHAIN'
18 include 'COMMON.INTERACT'
20 include 'COMMON.LOCAL'
21 include 'COMMON.TORSION'
22 include 'COMMON.SCCOR'
23 include 'COMMON.SCROT'
24 include 'COMMON.FFIELD'
25 include 'COMMON.NAMES'
26 include 'COMMON.SBRIDGE'
29 include 'COMMON.LAGRANGE.5diag'
31 include 'COMMON.LAGRANGE'
35 include 'COMMON.LANGEVIN.lang0.5diag'
37 include 'COMMON.LANGEVIN.lang0'
40 include 'COMMON.LANGEVIN'
42 include 'COMMON.SETUP'
43 include 'COMMON.CONTROL'
44 include 'COMMON.SHIELD'
46 character*1 onelett(4) /"G","A","P","D"/
47 character*1 toronelet(-2:2) /"p","a","G","A","P"/
49 double precision blower(3,3,maxlob)
51 character*3 lancuch,ucase
52 character*1000 weightcard
54 integer i,ii,j,jj,k,kk,l,ll,lll,llll,m,mm,n,iblock,junk,ijunk,
56 double precision akl,v0ij,si,rri,epsij,v0ijsccor,epsijlip,rjunk,
57 & sigt2sq,sigt1sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,
59 double precision dwa16
61 C For printing parameters after they are read set the following in the UNRES
64 C setenv PRINT_PARM YES
66 C To print parameters in LaTeX format rather than as ASCII tables:
70 call getenv_loc("PRINT_PARM",lancuch)
71 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
72 & .and. (me.eq.king.or..not.out1file) .and. fg_rank.eq.0
73 call getenv_loc("LATEX",lancuch)
74 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
76 dwa16=2.0d0**(1.0d0/6.0d0)
79 C Assign virtual-bond length
84 c Read the virtual-bond parameters, masses, and moments of inertia
85 c and Stokes radii of the peptide group and side chains
88 read (ibond,*) vbldp0,vbldpdum,akp,mp,ip,pstok
91 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
96 dsc_inv(i)=1.0D0/dsc(i)
100 read (ibond,*) junk,vbldp0,vbldpdum,akp,rjunk,mp,ip,pstok
102 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
103 & j=1,nbondterm(i)),msc(i),isc(i),restok(i)
104 dsc(i) = vbldsc0(1,i)
108 dsc_inv(i)=1.0D0/dsc(i)
113 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
114 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',
116 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
118 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
119 & vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
121 write (iout,'(13x,3f10.5)')
122 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
127 C reading lipid parameters
129 write (iout,*) "iliptranpar",iliptranpar
132 read(iliptranpar,*) pepliptran
134 read(iliptranpar,*) liptranene(i)
138 write (iout,'(/a)') "Water-lipid transfer parameters"
139 write (iout,'(a3,3x,f10.5)') 'p',pepliptran
141 write (iout,'(a3,3x,f10.5)') restyp(i),liptranene(i)
146 C Read the parameters of the probability distribution/energy expression
147 C of the virtual-bond valence angles theta
150 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),
151 & (bthet(j,i,1,1),j=1,2)
152 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
153 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
154 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
158 athet(1,i,1,-1)=athet(1,i,1,1)
159 athet(2,i,1,-1)=athet(2,i,1,1)
160 bthet(1,i,1,-1)=-bthet(1,i,1,1)
161 bthet(2,i,1,-1)=-bthet(2,i,1,1)
162 athet(1,i,-1,1)=-athet(1,i,1,1)
163 athet(2,i,-1,1)=-athet(2,i,1,1)
164 bthet(1,i,-1,1)=bthet(1,i,1,1)
165 bthet(2,i,-1,1)=bthet(2,i,1,1)
169 athet(1,i,-1,-1)=athet(1,-i,1,1)
170 athet(2,i,-1,-1)=-athet(2,-i,1,1)
171 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
172 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
173 athet(1,i,-1,1)=athet(1,-i,1,1)
174 athet(2,i,-1,1)=-athet(2,-i,1,1)
175 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
176 bthet(2,i,-1,1)=bthet(2,-i,1,1)
177 athet(1,i,1,-1)=-athet(1,-i,1,1)
178 athet(2,i,1,-1)=athet(2,-i,1,1)
179 bthet(1,i,1,-1)=bthet(1,-i,1,1)
180 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
185 polthet(j,i)=polthet(j,-i)
188 gthet(j,i)=gthet(j,-i)
196 & 'Parameters of the virtual-bond valence angles:'
197 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
198 & ' ATHETA0 ',' A1 ',' A2 ',
201 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
202 & a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
204 write (iout,'(/a/9x,5a/79(1h-))')
205 & 'Parameters of the expression for sigma(theta_c):',
206 & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
207 & ' ALPH3 ',' SIGMA0C '
209 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
210 & (polthet(j,i),j=0,3),sigc0(i)
212 write (iout,'(/a/9x,5a/79(1h-))')
213 & 'Parameters of the second gaussian:',
214 & ' THETA0 ',' SIGMA0 ',' G1 ',
217 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
218 & sig0(i),(gthet(j,i),j=1,3)
222 & 'Parameters of the virtual-bond valence angles:'
223 write (iout,'(/a/9x,5a/79(1h-))')
224 & 'Coefficients of expansion',
225 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
226 & ' b1*10^1 ',' b2*10^1 '
228 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
229 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
230 & (10*bthet(j,i,1,1),j=1,2)
232 write (iout,'(/a/9x,5a/79(1h-))')
233 & 'Parameters of the expression for sigma(theta_c):',
234 & ' alpha0 ',' alph1 ',' alph2 ',
235 & ' alhp3 ',' sigma0c '
237 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
238 & (polthet(j,i),j=0,3),sigc0(i)
240 write (iout,'(/a/9x,5a/79(1h-))')
241 & 'Parameters of the second gaussian:',
242 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
245 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
246 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
251 IF (tor_mode.eq.0) THEN
253 C Read the parameters of Utheta determined from ab initio surfaces
254 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
256 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,
257 & ntheterm3,nsingle,ndouble
258 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
259 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
261 ithetyp(i)=-ithetyp(-i)
264 do i=-maxthetyp,maxthetyp
265 do j=-maxthetyp,maxthetyp
266 do k=-maxthetyp,maxthetyp
267 aa0thet(i,j,k,iblock)=0.0d0
269 aathet(l,i,j,k,iblock)=0.0d0
273 bbthet(m,l,i,j,k,iblock)=0.0d0
274 ccthet(m,l,i,j,k,iblock)=0.0d0
275 ddthet(m,l,i,j,k,iblock)=0.0d0
276 eethet(m,l,i,j,k,iblock)=0.0d0
282 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
283 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
291 c VAR:iblock means terminally blocking group 1=non-proline 2=proline
293 c VAR:ntethtyp is type of theta potentials type currently 0=glycine
294 c VAR:1=non-glicyne non-proline 2=proline
295 c VAR:negative values for D-aminoacid
297 do j=-nthetyp,nthetyp
298 do k=-nthetyp,nthetyp
299 read (ithep,'(6a)',end=111,err=111) res1
300 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
301 c VAR: aa0thet is variable describing the average value of Foureir
302 c VAR: expansion series
303 c VAR: aathet is foureir expansion in theta/2 angle for full formula
304 c VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
305 Condens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
306 read (ithep,*,end=111,err=111)
307 &(aathet(l,i,j,k,iblock),l=1,ntheterm)
308 read (ithep,*,end=111,err=111)
309 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
310 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
311 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
312 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),
314 read (ithep,*,end=111,err=111)
315 & (((ffthet(llll,lll,ll,i,j,k,iblock),
316 & ffthet(lll,llll,ll,i,j,k,iblock),
317 & ggthet(llll,lll,ll,i,j,k,iblock),
318 & ggthet(lll,llll,ll,i,j,k,iblock),
319 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
324 C For dummy ends assign glycine-type coefficients of theta-only terms; the
325 C coefficients of theta-and-gamma-dependent terms are zero.
326 C IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
327 C RECOMENTDED AFTER VERSION 3.3)
331 c aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
332 c aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
334 c aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
335 c aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
338 c aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
340 c aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
343 C AND COMMENT THE LOOPS BELOW
347 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
348 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
350 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
351 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
354 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
356 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
360 C Substitution for D aminoacids from symmetry.
363 do j=-nthetyp,nthetyp
364 do k=-nthetyp,nthetyp
365 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
367 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
371 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
372 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
373 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
374 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
380 ffthet(llll,lll,ll,i,j,k,iblock)=
381 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
382 ffthet(lll,llll,ll,i,j,k,iblock)=
383 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
384 ggthet(llll,lll,ll,i,j,k,iblock)=
385 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
386 ggthet(lll,llll,ll,i,j,k,iblock)=
387 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
396 C Control printout of the coefficients of virtual-bond-angle potentials
399 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
402 do j=-nthetyp,nthetyp
403 do k=-nthetyp,nthetyp
404 write (iout,'(//4a)')
405 & 'Type ',toronelet(i),toronelet(j),toronelet(k)
406 write (iout,'(//a,10x,a)') " l","a[l]"
407 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
408 write (iout,'(i2,1pe15.5)')
409 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
411 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))')
412 & "b",l,"c",l,"d",l,"e",l
414 write (iout,'(i2,4(1pe15.5))') m,
415 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
416 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
420 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
421 & "f+",l,"f-",l,"g+",l,"g-",l
424 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
425 & ffthet(n,m,l,i,j,k,iblock),
426 & ffthet(m,n,l,i,j,k,iblock),
427 & ggthet(n,m,l,i,j,k,iblock),
428 & ggthet(m,n,l,i,j,k,iblock)
441 C here will be the apropriate recalibrating for D-aminoacid
442 read (ithep,*,end=121,err=121) nthetyp
443 do i=-nthetyp+1,nthetyp-1
444 read (ithep,*,end=121,err=121) nbend_kcc_Tb(i)
445 do j=0,nbend_kcc_Tb(i)
446 read (ithep,*,end=121,err=121) ijunk,v1bend_chyb(j,i)
451 & "Parameters of the valence-only potentials"
452 do i=-nthetyp+1,nthetyp-1
453 write (iout,'(2a)') "Type ",toronelet(i)
454 do k=0,nbend_kcc_Tb(i)
455 write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
462 c write (2,*) "Start reading THETA_PDB",ithep_pdb
465 read (ithep_pdb,*,err=111,end=111)
466 & a0thet(i),(athet(j,i,1,1),j=1,2),
467 & (bthet(j,i,1,1),j=1,2)
468 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
469 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
470 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
474 athet(1,i,1,-1)=athet(1,i,1,1)
475 athet(2,i,1,-1)=athet(2,i,1,1)
476 bthet(1,i,1,-1)=-bthet(1,i,1,1)
477 bthet(2,i,1,-1)=-bthet(2,i,1,1)
478 athet(1,i,-1,1)=-athet(1,i,1,1)
479 athet(2,i,-1,1)=-athet(2,i,1,1)
480 bthet(1,i,-1,1)=bthet(1,i,1,1)
481 bthet(2,i,-1,1)=bthet(2,i,1,1)
485 athet(1,i,-1,-1)=athet(1,-i,1,1)
486 athet(2,i,-1,-1)=-athet(2,-i,1,1)
487 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
488 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
489 athet(1,i,-1,1)=athet(1,-i,1,1)
490 athet(2,i,-1,1)=-athet(2,-i,1,1)
491 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
492 bthet(2,i,-1,1)=bthet(2,-i,1,1)
493 athet(1,i,1,-1)=-athet(1,-i,1,1)
494 athet(2,i,1,-1)=athet(2,-i,1,1)
495 bthet(1,i,1,-1)=bthet(1,-i,1,1)
496 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
501 polthet(j,i)=polthet(j,-i)
504 gthet(j,i)=gthet(j,-i)
507 c write (2,*) "End reading THETA_PDB"
513 C Read the parameters of the probability distribution/energy expression
514 C of the side chains.
517 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
521 dsc_inv(i)=1.0D0/dsc(i)
532 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
533 & ((blower(k,l,1),l=1,k),k=1,3)
534 censc(1,1,-i)=censc(1,1,i)
535 censc(2,1,-i)=censc(2,1,i)
536 censc(3,1,-i)=-censc(3,1,i)
538 read (irotam,*,end=112,err=112) bsc(j,i)
539 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
540 & ((blower(k,l,j),l=1,k),k=1,3)
541 censc(1,j,-i)=censc(1,j,i)
542 censc(2,j,-i)=censc(2,j,i)
543 censc(3,j,-i)=-censc(3,j,i)
544 C BSC is amplitude of Gaussian
551 akl=akl+blower(k,m,j)*blower(l,m,j)
555 if (((k.eq.3).and.(l.ne.3))
556 & .or.((l.eq.3).and.(k.ne.3))) then
557 gaussc(k,l,j,-i)=-akl
558 gaussc(l,k,j,-i)=-akl
570 write (iout,'(/a)') 'Parameters of side-chain local geometry'
575 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
576 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
577 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
578 & 'log h',(bsc(j,i),j=1,nlobi)
579 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
580 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
582 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
583 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
586 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
587 write (iout,'(a,f10.4,4(16x,f10.4))')
588 & 'Center ',(bsc(j,i),j=1,nlobi)
589 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
598 C Read scrot parameters for potentials determined from all-atom AM1 calculations
599 C added by Urszula Kozlowska 07/11/2007
602 read (irotam,*,end=112,err=112)
604 read (irotam,*,end=112,err=112)
607 read(irotam,*,end=112,err=112) sc_parmin(j,i)
612 C Read the parameters of the probability distribution/energy expression
613 C of the side chains.
615 write (2,*) "Start reading ROTAM_PDB"
617 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
621 dsc_inv(i)=1.0D0/dsc(i)
632 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
633 & ((blower(k,l,1),l=1,k),k=1,3)
635 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
636 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
637 & ((blower(k,l,j),l=1,k),k=1,3)
644 akl=akl+blower(k,m,j)*blower(l,m,j)
654 write (2,*) "End reading ROTAM_PDB"
658 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
659 C interaction energy of the Gly, Ala, and Pro prototypes.
661 read (ifourier,*) nloctyp
662 SPLIT_FOURIERTOR = nloctyp.lt.0
663 nloctyp = iabs(nloctyp)
665 read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
666 read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
667 itype2loc(ntyp1)=nloctyp
668 iloctyp(nloctyp)=ntyp1
670 itype2loc(-i)=-itype2loc(i)
679 iloctyp(-i)=-iloctyp(i)
681 c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
682 c write (iout,*) "nloctyp",nloctyp,
683 c & " iloctyp",(iloctyp(i),i=0,nloctyp)
692 c write (iout,*) "NEWCORR",i
693 read (ifourier,*,end=115,err=115)
696 read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
699 c write (iout,*) "NEWCORR BNEW1"
700 c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
703 read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
706 c write (iout,*) "NEWCORR BNEW2"
707 c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
709 read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
710 read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
712 c write (iout,*) "NEWCORR CCNEW"
713 c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
715 read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
716 read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
718 c write (iout,*) "NEWCORR DDNEW"
719 c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
723 read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
727 c write (iout,*) "NEWCORR EENEW1"
728 c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
730 read (ifourier,*,end=115,err=115) e0new(ii,i)
732 c write (iout,*) (e0new(ii,i),ii=1,3)
734 c write (iout,*) "NEWCORR EENEW"
737 ccnew(ii,1,i)=ccnew(ii,1,i)/2
738 ccnew(ii,2,i)=ccnew(ii,2,i)/2
739 ddnew(ii,1,i)=ddnew(ii,1,i)/2
740 ddnew(ii,2,i)=ddnew(ii,2,i)/2
745 bnew1(ii,1,-i)= bnew1(ii,1,i)
746 bnew1(ii,2,-i)=-bnew1(ii,2,i)
747 bnew2(ii,1,-i)= bnew2(ii,1,i)
748 bnew2(ii,2,-i)=-bnew2(ii,2,i)
751 c ccnew(ii,1,i)=ccnew(ii,1,i)/2
752 c ccnew(ii,2,i)=ccnew(ii,2,i)/2
753 c ddnew(ii,1,i)=ddnew(ii,1,i)/2
754 c ddnew(ii,2,i)=ddnew(ii,2,i)/2
755 ccnew(ii,1,-i)=ccnew(ii,1,i)
756 ccnew(ii,2,-i)=-ccnew(ii,2,i)
757 ddnew(ii,1,-i)=ddnew(ii,1,i)
758 ddnew(ii,2,-i)=-ddnew(ii,2,i)
760 e0new(1,-i)= e0new(1,i)
761 e0new(2,-i)=-e0new(2,i)
762 e0new(3,-i)=-e0new(3,i)
764 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
765 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
766 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
767 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
771 write (iout,'(a)') "Coefficients of the multibody terms"
772 c do i=-nloctyp+1,nloctyp-1
773 do i=-nloctyp,nloctyp
774 write (iout,*) "Type: ",onelet(iloctyp(i))
775 write (iout,*) "Coefficients of the expansion of B1"
777 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
779 write (iout,*) "Coefficients of the expansion of B2"
781 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
783 write (iout,*) "Coefficients of the expansion of C"
784 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
785 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
786 write (iout,*) "Coefficients of the expansion of D"
787 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
788 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
789 write (iout,*) "Coefficients of the expansion of E"
790 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
793 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
798 IF (SPLIT_FOURIERTOR) THEN
800 c write (iout,*) "NEWCORR TOR",i
801 read (ifourier,*,end=115,err=115)
804 read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i)
807 c write (iout,*) "NEWCORR BNEW1 TOR"
808 c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
811 read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i)
814 c write (iout,*) "NEWCORR BNEW2 TOR"
815 c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
817 read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i)
818 read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i)
820 c write (iout,*) "NEWCORR CCNEW TOR"
821 c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
823 read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i)
824 read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i)
826 c write (iout,*) "NEWCORR DDNEW TOR"
827 c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
831 read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i)
835 c write (iout,*) "NEWCORR EENEW1 TOR"
836 c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
838 read (ifourier,*,end=115,err=115) e0newtor(ii,i)
840 c write (iout,*) (e0newtor(ii,i),ii=1,3)
842 c write (iout,*) "NEWCORR EENEW TOR"
845 ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
846 ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
847 ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
848 ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
853 bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
854 bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
855 bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
856 bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
859 ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
860 ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
861 ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
862 ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
864 e0newtor(1,-i)= e0newtor(1,i)
865 e0newtor(2,-i)=-e0newtor(2,i)
866 e0newtor(3,-i)=-e0newtor(3,i)
868 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
869 eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
870 eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
871 eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
876 & "Single-body coefficients of the torsional potentials"
877 do i=-nloctyp+1,nloctyp-1
878 write (iout,*) "Type: ",onelet(iloctyp(i))
879 write (iout,*) "Coefficients of the expansion of B1tor"
881 write (iout,'(3hB1(,i1,1h),3f10.5)')
882 & j,(bnew1tor(k,j,i),k=1,3)
884 write (iout,*) "Coefficients of the expansion of B2tor"
886 write (iout,'(3hB2(,i1,1h),3f10.5)')
887 & j,(bnew2tor(k,j,i),k=1,3)
889 write (iout,*) "Coefficients of the expansion of Ctor"
890 write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
891 write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
892 write (iout,*) "Coefficients of the expansion of Dtor"
893 write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
894 write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
895 write (iout,*) "Coefficients of the expansion of Etor"
896 write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
899 write (iout,'(1hE,2i1,2f10.5)')
900 & j,k,(eenewtor(l,j,k,i),l=1,2)
906 do i=-nloctyp+1,nloctyp-1
909 bnew1tor(ii,j,i)=bnew1(ii,j,i)
914 bnew2tor(ii,j,i)=bnew2(ii,j,i)
918 ccnewtor(ii,1,i)=ccnew(ii,1,i)
919 ccnewtor(ii,2,i)=ccnew(ii,2,i)
920 ddnewtor(ii,1,i)=ddnew(ii,1,i)
921 ddnewtor(ii,2,i)=ddnew(ii,2,i)
924 ENDIF !SPLIT_FOURIER_TOR
927 & write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
929 read (ifourier,*,end=115,err=115)
930 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
932 write (iout,*) 'Type ',onelet(iloctyp(i))
933 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
947 c B1tilde(1,i) = b(3)
948 c B1tilde(2,i) =-b(5)
949 c B1tilde(1,-i) =-b(3)
950 c B1tilde(2,-i) =b(5)
957 cc B1tilde(1,i) = b(3,i)
958 cc B1tilde(2,i) =-b(5,i)
959 C B1tilde(1,-i) =-b(3,i)
960 C B1tilde(2,-i) =b(5,i)
961 cc b1tilde(1,i)=0.0d0
962 cc b1tilde(2,i)=0.0d0
974 CCold(1,1,-i)= b(7,i)
975 CCold(2,2,-i)=-b(7,i)
976 CCold(2,1,-i)=-b(9,i)
977 CCold(1,2,-i)=-b(9,i)
982 c Ctilde(1,1,i)= CCold(1,1,i)
983 c Ctilde(1,2,i)= CCold(1,2,i)
984 c Ctilde(2,1,i)=-CCold(2,1,i)
985 c Ctilde(2,2,i)=-CCold(2,2,i)
987 c Ctilde(1,1,i)=0.0d0
988 c Ctilde(1,2,i)=0.0d0
989 c Ctilde(2,1,i)=0.0d0
990 c Ctilde(2,2,i)=0.0d0
995 DDold(1,1,-i)= b(6,i)
996 DDold(2,2,-i)=-b(6,i)
997 DDold(2,1,-i)=-b(8,i)
998 DDold(1,2,-i)=-b(8,i)
1003 c Dtilde(1,1,i)= DD(1,1,i)
1004 c Dtilde(1,2,i)= DD(1,2,i)
1005 c Dtilde(2,1,i)=-DD(2,1,i)
1006 c Dtilde(2,2,i)=-DD(2,2,i)
1008 c Dtilde(1,1,i)=0.0d0
1009 c Dtilde(1,2,i)=0.0d0
1010 c Dtilde(2,1,i)=0.0d0
1011 c Dtilde(2,2,i)=0.0d0
1012 EEold(1,1,i)= b(10,i)+b(11,i)
1013 EEold(2,2,i)=-b(10,i)+b(11,i)
1014 EEold(2,1,i)= b(12,i)-b(13,i)
1015 EEold(1,2,i)= b(12,i)+b(13,i)
1016 EEold(1,1,-i)= b(10,i)+b(11,i)
1017 EEold(2,2,-i)=-b(10,i)+b(11,i)
1018 EEold(2,1,-i)=-b(12,i)+b(13,i)
1019 EEold(1,2,-i)=-b(12,i)-b(13,i)
1020 write(iout,*) "TU DOCHODZE"
1026 c ee(2,1,i)=ee(1,2,i)
1031 &"Coefficients of the cumulants (independent of valence angles)"
1032 do i=-nloctyp+1,nloctyp-1
1033 write (iout,*) 'Type ',onelet(iloctyp(i))
1035 write(iout,'(2f10.5)') B(3,i),B(5,i)
1037 write(iout,'(2f10.5)') B(2,i),B(4,i)
1040 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
1044 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
1048 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1056 C Read torsional parameters in old format
1058 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1059 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1060 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1065 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1071 write (iout,'(/a/)') 'Torsional constants:'
1074 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1075 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1081 C Read torsional parameters
1083 IF (TOR_MODE.eq.0) THEN
1085 read (itorp,*,end=113,err=113) ntortyp
1086 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1088 itype2loc(i)=itortyp(i)
1091 itype2loc(-i)=-itype2loc(i)
1093 itortyp(ntyp1)=ntortyp
1096 itortyp(i)=-itortyp(-i)
1098 write (iout,*) 'ntortyp',ntortyp
1100 do j=-ntortyp+1,ntortyp-1
1101 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
1103 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1104 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1107 do k=1,nterm(i,j,iblock)
1108 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
1110 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1111 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1112 v0ij=v0ij+si*v1(k,i,j,iblock)
1114 c write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
1115 c write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
1116 c &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
1118 do k=1,nlor(i,j,iblock)
1119 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
1120 & vlor2(k,i,j),vlor3(k,i,j)
1121 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1124 v0(-i,-j,iblock)=v0ij
1130 write (iout,'(/a/)') 'Parameters of the SCCOR potentials:'
1133 do j=-ntortyp+1,ntortyp-1
1134 write (iout,*) 'ityp',i,' jtyp',j
1135 write (iout,*) 'Fourier constants'
1136 do k=1,nterm(i,j,iblock)
1137 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
1140 write (iout,*) 'Lorenz constants'
1141 do k=1,nlor(i,j,iblock)
1142 write (iout,'(3(1pe15.5))')
1143 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1151 C 6/23/01 Read parameters for double torsionals
1155 do j=-ntortyp+1,ntortyp-1
1156 do k=-ntortyp+1,ntortyp-1
1157 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1158 c write (iout,*) "OK onelett",
1161 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
1162 & .or. t3.ne.toronelet(k)) then
1163 write (iout,*) "Error in double torsional parameter file",
1166 call MPI_Finalize(Ierror)
1168 stop "Error in double torsional parameter file"
1170 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
1171 & ntermd_2(i,j,k,iblock)
1172 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1173 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1174 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
1175 & ntermd_1(i,j,k,iblock))
1176 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
1177 & ntermd_1(i,j,k,iblock))
1178 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
1179 & ntermd_1(i,j,k,iblock))
1180 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
1181 & ntermd_1(i,j,k,iblock))
1182 C Martix of D parameters for one dimesional foureir series
1183 do l=1,ntermd_1(i,j,k,iblock)
1184 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1185 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1186 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1187 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1188 c write(iout,*) "whcodze" ,
1189 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1191 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
1192 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
1193 & v2s(m,l,i,j,k,iblock),
1194 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1195 C Martix of D parameters for two dimesional fourier series
1196 do l=1,ntermd_2(i,j,k,iblock)
1198 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1199 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1200 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1201 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1210 write (iout,*) 'Constants for double torsionals'
1213 do j=-ntortyp+1,ntortyp-1
1214 do k=-ntortyp+1,ntortyp-1
1215 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
1216 & ' nsingle',ntermd_1(i,j,k,iblock),
1217 & ' ndouble',ntermd_2(i,j,k,iblock)
1219 write (iout,*) 'Single angles:'
1220 do l=1,ntermd_1(i,j,k,iblock)
1221 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
1222 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
1223 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
1224 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1227 write (iout,*) 'Pairs of angles:'
1228 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1229 do l=1,ntermd_2(i,j,k,iblock)
1230 write (iout,'(i5,20f10.5)')
1231 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1234 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1235 do l=1,ntermd_2(i,j,k,iblock)
1236 write (iout,'(i5,20f10.5)')
1237 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
1238 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1247 ELSE IF (TOR_MODE.eq.1) THEN
1249 C read valence-torsional parameters
1250 read (itorp,*,end=121,err=121) ntortyp
1252 write (iout,*) "Valence-torsional parameters read in ntortyp",
1254 read (itorp,*,end=121,err=121) (itortyp(i),i=1,ntyp)
1255 write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
1258 itype2loc(i)=itortyp(i)
1262 itortyp(i)=-itortyp(-i)
1264 do i=-ntortyp+1,ntortyp-1
1265 do j=-ntortyp+1,ntortyp-1
1266 C first we read the cos and sin gamma parameters
1267 read (itorp,'(13x,a)',end=121,err=121) string
1268 write (iout,*) i,j,string
1269 read (itorp,*,end=121,err=121)
1270 & nterm_kcc(j,i),nterm_kcc_Tb(j,i)
1271 C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
1272 do k=1,nterm_kcc(j,i)
1273 do l=1,nterm_kcc_Tb(j,i)
1274 do ll=1,nterm_kcc_Tb(j,i)
1275 read (itorp,*,end=121,err=121) ii,jj,kk,
1276 & v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1284 c AL 4/8/16: Calculate coefficients from one-body parameters
1287 itortyp(i)=itype2loc(i)
1290 &"Val-tor parameters calculated from cumulant coefficients ntortyp"
1292 do i=-ntortyp+1,ntortyp-1
1293 do j=-ntortyp+1,ntortyp-1
1296 do k=1,nterm_kcc_Tb(j,i)
1297 do l=1,nterm_kcc_Tb(j,i)
1298 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)
1299 & +bnew1tor(k,2,i)*bnew2tor(l,2,j)
1300 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)
1301 & +bnew1tor(k,2,i)*bnew2tor(l,1,j)
1304 do k=1,nterm_kcc_Tb(j,i)
1305 do l=1,nterm_kcc_Tb(j,i)
1307 v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1308 & -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1309 v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1310 & +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1312 v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1313 & -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1314 v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1315 & +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1319 c f(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)
1323 write (iout,*) "TOR_MODE>1 only with NEWCORR"
1328 if (tor_mode.gt.0 .and. lprint) then
1329 c Print valence-torsional parameters
1331 & "Parameters of the valence-torsional potentials"
1332 do i=-ntortyp+1,ntortyp-1
1333 do j=-ntortyp+1,ntortyp-1
1334 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
1335 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
1336 do k=1,nterm_kcc(j,i)
1337 do l=1,nterm_kcc_Tb(j,i)
1338 do ll=1,nterm_kcc_Tb(j,i)
1339 write (iout,'(3i5,2f15.4)')
1340 & k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1349 C Read of Side-chain backbone correlation parameters
1350 C Modified 11 May 2012 by Adasko
1353 read (isccor,*,end=119,err=119) nsccortyp
1355 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1357 isccortyp(i)=-isccortyp(-i)
1359 iscprol=isccortyp(20)
1360 c write (iout,*) 'ntortyp',ntortyp
1362 cc maxinter is maximum interaction sites
1366 read (isccor,*,end=119,err=119)
1367 &nterm_sccor(i,j),nlor_sccor(i,j)
1373 nterm_sccor(-i,j)=nterm_sccor(i,j)
1374 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1375 nterm_sccor(i,-j)=nterm_sccor(i,j)
1376 do k=1,nterm_sccor(i,j)
1377 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
1379 if (j.eq.iscprol) then
1380 if (i.eq.isccortyp(10)) then
1381 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1382 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1384 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
1385 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1386 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
1387 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1388 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1389 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1390 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1391 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1394 if (i.eq.isccortyp(10)) then
1395 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1396 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1398 if (j.eq.isccortyp(10)) then
1399 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1400 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1402 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1403 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1404 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1405 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1406 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1407 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1411 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1412 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1413 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1414 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1417 do k=1,nlor_sccor(i,j)
1418 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
1419 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1420 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
1421 &(1+vlor3sccor(k,i,j)**2)
1423 v0sccor(l,i,j)=v0ijsccor
1424 v0sccor(l,-i,j)=v0ijsccor1
1425 v0sccor(l,i,-j)=v0ijsccor2
1426 v0sccor(l,-i,-j)=v0ijsccor3
1432 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1433 c write (iout,*) 'ntortyp',ntortyp
1435 cc maxinter is maximum interaction sites
1439 read (isccor,*,end=119,err=119)
1440 & nterm_sccor(i,j),nlor_sccor(i,j)
1444 do k=1,nterm_sccor(i,j)
1445 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
1447 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1450 do k=1,nlor_sccor(i,j)
1451 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
1452 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1453 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
1454 &(1+vlor3sccor(k,i,j)**2)
1456 v0sccor(l,i,j)=v0ijsccor
1464 write (iout,'(/a/)') 'SCCor torsional constants:'
1468 write (iout,*) 'ityp',i,' jtyp',j
1469 write (iout,*) 'Fourier constants'
1470 do k=1,nterm_sccor(i,j)
1471 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1473 write (iout,*) 'Lorenz constants'
1474 do k=1,nlor_sccor(i,j)
1475 write (iout,'(3(1pe15.5))')
1476 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1485 C Read electrostatic-interaction parameters
1489 write (iout,'(/a)') 'Electrostatic interaction constants:'
1490 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1491 & 'IT','JT','APP','BPP','AEL6','AEL3'
1493 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1494 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1495 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1496 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1501 app (i,j)=epp(i,j)*rri*rri
1502 bpp (i,j)=-2.0D0*epp(i,j)*rri
1503 ael6(i,j)=elpp6(i,j)*4.2D0**6
1504 ael3(i,j)=elpp3(i,j)*4.2D0**3
1506 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1507 & ael6(i,j),ael3(i,j)
1512 C Read side-chain interaction parameters.
1514 read (isidep,*,end=117,err=117) ipot,expon
1515 if (ipot.lt.1 .or. ipot.gt.5) then
1516 write (iout,'(2a)') 'Error while reading SC interaction',
1517 & 'potential file - unknown potential type.'
1519 call MPI_Finalize(Ierror)
1525 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1526 & ', exponents are ',expon,2*expon
1527 goto (10,20,30,30,40) ipot
1528 C----------------------- LJ potential ---------------------------------
1529 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1530 & (sigma0(i),i=1,ntyp)
1532 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1533 write (iout,'(a/)') 'The epsilon array:'
1534 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1535 write (iout,'(/a)') 'One-body parameters:'
1536 write (iout,'(a,4x,a)') 'residue','sigma'
1537 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1540 C----------------------- LJK potential --------------------------------
1541 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1542 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1544 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1545 write (iout,'(a/)') 'The epsilon array:'
1546 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1547 write (iout,'(/a)') 'One-body parameters:'
1548 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1549 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1553 C---------------------- GB or BP potential -----------------------------
1555 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1557 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
1558 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
1559 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
1560 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
1561 C now we start reading lipid
1563 read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
1564 if (lprint) write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp)
1566 C print *,"WARNING!!"
1568 C epslip(i,j)=epslip(i,j)+0.05d0
1571 if (lprint) write(iout,*) epslip(1,1),"OK?"
1572 C For the GB potential convert sigma'**2 into chi'
1575 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1579 write (iout,'(/a/)') 'Parameters of the BP potential:'
1580 write (iout,'(a/)') 'The epsilon array:'
1581 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1582 write (iout,'(/a)') 'One-body parameters:'
1583 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1585 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1586 & chip(i),alp(i),i=1,ntyp)
1589 C--------------------- GBV potential -----------------------------------
1590 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1591 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1592 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1594 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1595 write (iout,'(a/)') 'The epsilon array:'
1596 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1597 write (iout,'(/a)') 'One-body parameters:'
1598 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1599 & 's||/s_|_^2',' chip ',' alph '
1600 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1601 & sigii(i),chip(i),alp(i),i=1,ntyp)
1605 C-----------------------------------------------------------------------
1606 C Calculate the "working" parameters of SC interactions.
1610 epslip(i,j)=epslip(j,i)
1615 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1616 sigma(j,i)=sigma(i,j)
1617 rs0(i,j)=dwa16*sigma(i,j)
1621 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1622 & 'Working parameters of the SC interactions:',
1623 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1628 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1637 sigeps=dsign(1.0D0,epsij)
1639 aa_aq(i,j)=epsij*rrij*rrij
1640 bb_aq(i,j)=-sigeps*epsij*rrij
1641 aa_aq(j,i)=aa_aq(i,j)
1642 bb_aq(j,i)=bb_aq(i,j)
1643 epsijlip=epslip(i,j)
1644 sigeps=dsign(1.0D0,epsijlip)
1645 epsijlip=dabs(epsijlip)
1646 aa_lip(i,j)=epsijlip*rrij*rrij
1647 bb_lip(i,j)=-sigeps*epsijlip*rrij
1648 aa_lip(j,i)=aa_lip(i,j)
1649 bb_lip(j,i)=bb_lip(i,j)
1651 sigt1sq=sigma0(i)**2
1652 sigt2sq=sigma0(j)**2
1655 ratsig1=sigt2sq/sigt1sq
1656 ratsig2=1.0D0/ratsig1
1657 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1658 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1659 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1663 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1664 sigmaii(i,j)=rsum_max
1665 sigmaii(j,i)=rsum_max
1667 c sigmaii(i,j)=r0(i,j)
1668 c sigmaii(j,i)=r0(i,j)
1670 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1671 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1672 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1673 augm(i,j)=epsij*r_augm**(2*expon)
1674 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1681 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1682 & restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1683 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1688 write(iout,*) "tube param"
1689 read(itube,*) epspeptube,sigmapeptube
1690 sigmapeptube=sigmapeptube**6
1691 sigeps=dsign(1.0D0,epspeptube)
1692 epspeptube=dabs(epspeptube)
1693 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
1694 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
1695 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
1697 read(itube,*) epssctube,sigmasctube
1698 sigmasctube=sigmasctube**6
1699 sigeps=dsign(1.0D0,epssctube)
1700 epssctube=dabs(epssctube)
1701 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
1702 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
1703 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
1709 C Define the SC-p interaction constants (hard-coded; old style)
1712 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1714 c aad(i,1)=0.3D0*4.0D0**12
1715 C Following line for constants currently implemented
1716 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1717 aad(i,1)=1.5D0*4.0D0**12
1718 c aad(i,1)=0.17D0*5.6D0**12
1720 C "Soft" SC-p repulsion
1722 C Following line for constants currently implemented
1723 c aad(i,1)=0.3D0*4.0D0**6
1724 C "Hard" SC-p repulsion
1725 bad(i,1)=3.0D0*4.0D0**6
1726 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1735 C 8/9/01 Read the SC-p interaction constants from file
1738 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1741 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1742 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1743 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1744 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1748 write (iout,*) "Parameters of SC-p interactions:"
1750 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1751 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1757 C Define the constants of the disulfide bridge
1761 c Old arbitrary potential - commented out.
1766 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1767 c energy surface of diethyl disulfide.
1768 c A. Liwo and U. Kozlowska, 11/24/03
1784 C if(me.eq.king) then
1785 C write (iout,'(/a)') "Disulfide bridge parameters:"
1786 C write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1787 C write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1788 C write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1789 C write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1792 if (shield_mode.gt.0) then
1793 C VSolvSphere the volume of solving sphere
1794 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1795 C there will be no distinction between proline peptide group and normal peptide
1796 C group in case of shielding parameters
1797 c write (iout,*) "rpp(1,1)",rpp(1,1)," pi",pi
1798 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1799 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1800 write (iout,*) "VSolvSphere",VSolvSphere,"VSolvSphere_div",
1802 C long axis of side chain
1804 long_r_sidechain(i)=vbldsc0(1,i)
1805 short_r_sidechain(i)=sigma0(i)
1810 111 write (iout,*) "Error reading bending energy parameters."
1812 112 write (iout,*) "Error reading rotamer energy parameters."
1814 113 write (iout,*) "Error reading torsional energy parameters."
1816 114 write (iout,*) "Error reading double torsional energy parameters."
1819 & "Error reading cumulant (multibody energy) parameters."
1821 116 write (iout,*) "Error reading electrostatic energy parameters."
1823 1161 write (iout,*) "Error reading electrostatic energy parameters.Lip"
1825 117 write (iout,*) "Error reading side chain interaction parameters."
1827 118 write (iout,*) "Error reading SCp interaction parameters."
1829 119 write (iout,*) "Error reading SCCOR parameters"
1831 121 write (iout,*) "Error in Czybyshev parameters"
1834 call MPI_Finalize(Ierror)
1841 subroutine getenv_loc(var, val)
1842 character(*) var, val
1845 character(2000) line
1848 open (196,file='env',status='old',readonly,shared)
1850 c write(*,*)'looking for ',var
1851 10 read(196,*,err=11,end=11)line
1852 iread=index(line,var)
1853 c write(*,*)iread,' ',var,' ',line
1854 if (iread.eq.0) go to 10
1855 c write(*,*)'---> ',line
1861 iread=iread+ilen(var)+1
1862 read (line(iread:),*,err=12,end=12) val
1863 c write(*,*)'OK: ',var,' = ',val
1869 #elif (defined CRAY)
1870 integer lennam,lenval,ierror
1872 c getenv using a POSIX call, useful on the T3D
1873 c Sept 1996, comment out error check on advice of H. Pritchard
1876 if(lennam.le.0) stop '--error calling getenv--'
1877 call pxfgetenv(var,lennam,val,lenval,ierror)
1878 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1880 call getenv(var,val)
1882 C set the variables used for shielding effect
1883 C if (shield_mode.gt.0) then
1884 C VSolvSphere the volume of solving sphere
1886 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1887 C there will be no distinction between proline peptide group and normal peptide
1888 C group in case of shielding parameters
1889 C VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1890 C VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1891 C long axis of side chain
1893 C long_r_sidechain(i)=vbldsc0(1,i)
1894 C short_r_sidechain(i)=sigma0(i)
1896 C lets set the buffor value
1901 c---------------------------------------------------------------------------
1902 subroutine weightread
1904 C Read the parameters of the probability distributions of the virtual-bond
1905 C valence angles and the side chains and energy parameters.
1907 C Important! Energy-term weights ARE NOT read here; they are read from the
1908 C main input file instead, because NO defaults have yet been set for these
1912 include 'DIMENSIONS'
1917 include 'COMMON.IOUNITS'
1918 include 'COMMON.CHAIN'
1919 include 'COMMON.INTERACT'
1920 include 'COMMON.GEO'
1921 include 'COMMON.LOCAL'
1922 include 'COMMON.TORSION'
1923 include 'COMMON.SCCOR'
1924 include 'COMMON.SCROT'
1925 include 'COMMON.FFIELD'
1926 include 'COMMON.NAMES'
1927 include 'COMMON.SBRIDGE'
1929 include 'COMMON.SETUP'
1930 include 'COMMON.CONTROL'
1931 include 'COMMON.SHIELD'
1932 character*1000 weightcard
1934 double precision scalscp,wlong
1936 c READ energy-term weights
1938 call card_concat(weightcard)
1939 call reada(weightcard,'WLONG',wlong,1.0D0)
1940 call reada(weightcard,'WSC',wsc,wlong)
1941 call reada(weightcard,'WSCP',wscp,wlong)
1942 call reada(weightcard,'WELEC',welec,1.0D0)
1943 call reada(weightcard,'WVDWPP',wvdwpp,welec)
1944 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
1945 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
1946 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
1947 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
1948 call reada(weightcard,'WTURN3',wturn3,1.0D0)
1949 call reada(weightcard,'WTURN4',wturn4,1.0D0)
1950 call reada(weightcard,'WTURN6',wturn6,1.0D0)
1951 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
1952 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
1953 call reada(weightcard,'WBOND',wbond,1.0D0)
1954 call reada(weightcard,'WTOR',wtor,1.0D0)
1955 call reada(weightcard,'WTORD',wtor_d,1.0D0)
1956 call reada(weightcard,'WANG',wang,1.0D0)
1957 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
1958 call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
1959 call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
1960 call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
1961 call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
1962 call reada(weightcard,'SCAL14',scal14,0.4D0)
1963 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
1964 call reada(weightcard,'LIPSCALE',lipscale,1.0D0)
1965 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
1966 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
1967 call reada(weightcard,'TEMP0',temp0,300.0d0)
1968 call reada(weightcard,'WSHIELD',wshield,1.0d0)
1969 call reada(weightcard,'WLT',wliptran,0.0D0)
1970 call reada(weightcard,'WTUBE',wtube,1.0D0)
1971 call reada(weightcard,'WSAXS',wsaxs,1.0D0)
1972 if (index(weightcard,'SOFT').gt.0) ipot=6
1973 C 12/1/95 Added weight for the multi-body term WCORR
1974 call reada(weightcard,'WCORRH',wcorr,1.0D0)
1975 if (wcorr4.gt.0.0d0) wcorr=wcorr4
1996 weights(28)=wdfa_dist
1997 weights(29)=wdfa_tor
1998 weights(30)=wdfa_nei
1999 weights(31)=wdfa_beta
2001 aad(i,1)=scalscp*aad(i,1)
2002 aad(i,2)=scalscp*aad(i,2)
2003 bad(i,1)=scalscp*bad(i,1)
2004 bad(i,2)=scalscp*bad(i,2)
2006 if(me.eq.king.or..not.out1file)
2007 & write (iout,100) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
2008 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
2010 100 format (/'Energy-term weights (unscaled):'//
2011 & 'WSCC= ',f10.6,' (SC-SC)'/
2012 & 'WSCP= ',f10.6,' (SC-p)'/
2013 & 'WELEC= ',f10.6,' (p-p electr)'/
2014 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
2015 & 'WBOND= ',f10.6,' (stretching)'/
2016 & 'WANG= ',f10.6,' (bending)'/
2017 & 'WSCLOC= ',f10.6,' (SC local)'/
2018 & 'WTOR= ',f10.6,' (torsional)'/
2019 & 'WTORD= ',f10.6,' (double torsional)'/
2020 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
2021 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
2022 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
2023 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
2024 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
2025 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
2026 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
2027 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
2028 & 'WTURN6= ',f10.6,' (turns, 6th order)')
2029 if(me.eq.king.or..not.out1file)then
2030 if (wcorr4.gt.0.0d0) then
2031 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
2032 & 'between contact pairs of peptide groups'
2033 write (iout,'(2(a,f5.3/))')
2034 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
2035 & 'Range of quenching the correlation terms:',2*delt_corr
2036 else if (wcorr.gt.0.0d0) then
2037 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
2038 & 'between contact pairs of peptide groups'
2040 write (iout,'(a,f8.3)')
2041 & 'Scaling factor of 1,4 SC-p interactions:',scal14
2042 write (iout,'(a,f8.3)')
2043 & 'General scaling factor of SC-p interactions:',scalscp
2045 r0_corr=cutoff_corr-delt_corr
2047 call rescale_weights(t_bath)
2048 if(me.eq.king.or..not.out1file)
2049 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
2050 & wtor_d,wstrain,wel_loc,
2052 & wcorr,wcorr5,wcorr6,
2059 22 format (/'Energy-term weights (scaled):'//
2060 & 'WSCC= ',f10.6,' (SC-SC)'/
2061 & 'WSCP= ',f10.6,' (SC-p)'/
2062 & 'WELEC= ',f10.6,' (p-p electr)'/
2063 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
2064 & 'WBOND= ',f10.6,' (stretching)'/
2065 & 'WANG= ',f10.6,' (bending)'/
2066 & 'WSCLOC= ',f10.6,' (SC local)'/
2067 & 'WTOR= ',f10.6,' (torsional)'/
2068 & 'WTORD= ',f10.6,' (double torsional)'/
2069 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
2070 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
2072 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
2073 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
2074 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
2076 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
2077 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
2078 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
2080 & 'WTURN6= ',f10.6,' (turns, 6th order)'
2083 if(me.eq.king.or..not.out1file)
2084 & write (iout,*) "Reference temperature for weights calculation:",
2087 write (iout,'(/a/)') "DFA pseudopotential parameters:"
2088 write (iout,'(a,f10.6,a)')
2089 & "WDFAD= ",wdfa_dist," (distance)",
2090 & "WDFAT= ",wdfa_tor," (backbone angles)",
2091 & "WDFAN= ",wdfa_nei," (neighbors)",
2092 & "WDFAB= ",wdfa_beta," (beta structure)"
2094 call reada(weightcard,"D0CM",d0cm,3.78d0)
2095 call reada(weightcard,"AKCM",akcm,15.1d0)
2096 call reada(weightcard,"AKTH",akth,11.0d0)
2097 call reada(weightcard,"AKCT",akct,12.0d0)
2098 call reada(weightcard,"V1SS",v1ss,-1.08d0)
2099 call reada(weightcard,"V2SS",v2ss,7.61d0)
2100 call reada(weightcard,"V3SS",v3ss,13.7d0)
2101 call reada(weightcard,"EBR",ebr,-5.50D0)
2102 call reada(weightcard,"ATRISS",atriss,0.301D0)
2103 call reada(weightcard,"BTRISS",btriss,0.021D0)
2104 call reada(weightcard,"CTRISS",ctriss,1.001D0)
2105 call reada(weightcard,"DTRISS",dtriss,1.001D0)
2106 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
2108 dyn_ss_mask(i)=.false.
2112 dyn_ssbond_ij(i,j)=1.0d300
2116 call reada(weightcard,"HT",Ht,0.0D0)
2118 ss_depth=ebr/wsc-0.25*eps(1,1)
2119 Ht=Ht/wsc-0.25*eps(1,1)
2120 akcm=akcm*wstrain/wsc
2121 akth=akth*wstrain/wsc
2122 akct=akct*wstrain/wsc
2123 v1ss=v1ss*wstrain/wsc
2124 v2ss=v2ss*wstrain/wsc
2125 v3ss=v3ss*wstrain/wsc
2127 if (wstrain.ne.0.0) then
2128 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
2133 write (iout,*) "wshield,", wshield
2134 if(me.eq.king.or..not.out1file) then
2135 write (iout,*) "Parameters of the SS-bond potential:"
2136 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
2138 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
2139 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
2140 write (iout,*)" HT",Ht
2141 write (iout,*) "Parameters of the 'trisulfide' potential"
2142 write (iout,*) "ATRISS=", atriss
2143 write (iout,*) "BTRISS=", btriss
2144 write (iout,*) "CTRISS=", ctriss
2145 write (iout,*) "DTRISS=", dtriss
2146 c print *,'indpdb=',indpdb,' pdbref=',pdbref