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 call getenv_loc("LATEX",lancuch)
73 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
75 dwa16=2.0d0**(1.0d0/6.0d0)
78 C Assign virtual-bond length
83 c Read the virtual-bond parameters, masses, and moments of inertia
84 c and Stokes radii of the peptide group and side chains
87 read (ibond,*) vbldp0,vbldpdum,akp,mp,ip,pstok
90 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
95 dsc_inv(i)=1.0D0/dsc(i)
99 read (ibond,*) junk,vbldp0,vbldpdum,akp,rjunk,mp,ip,pstok
101 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
102 & j=1,nbondterm(i)),msc(i),isc(i),restok(i)
103 dsc(i) = vbldsc0(1,i)
107 dsc_inv(i)=1.0D0/dsc(i)
112 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
113 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',
115 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
117 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
118 & vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
120 write (iout,'(13x,3f10.5)')
121 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
125 C reading lipid parameters
127 write (iout,*) "iliptranpar",iliptranpar
130 read(iliptranpar,*) pepliptran
132 read(iliptranpar,*) liptranene(i)
137 C Read the parameters of the probability distribution/energy expression
138 C of the virtual-bond valence angles theta
141 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),
142 & (bthet(j,i,1,1),j=1,2)
143 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
144 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
145 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
149 athet(1,i,1,-1)=athet(1,i,1,1)
150 athet(2,i,1,-1)=athet(2,i,1,1)
151 bthet(1,i,1,-1)=-bthet(1,i,1,1)
152 bthet(2,i,1,-1)=-bthet(2,i,1,1)
153 athet(1,i,-1,1)=-athet(1,i,1,1)
154 athet(2,i,-1,1)=-athet(2,i,1,1)
155 bthet(1,i,-1,1)=bthet(1,i,1,1)
156 bthet(2,i,-1,1)=bthet(2,i,1,1)
160 athet(1,i,-1,-1)=athet(1,-i,1,1)
161 athet(2,i,-1,-1)=-athet(2,-i,1,1)
162 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
163 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
164 athet(1,i,-1,1)=athet(1,-i,1,1)
165 athet(2,i,-1,1)=-athet(2,-i,1,1)
166 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
167 bthet(2,i,-1,1)=bthet(2,-i,1,1)
168 athet(1,i,1,-1)=-athet(1,-i,1,1)
169 athet(2,i,1,-1)=athet(2,-i,1,1)
170 bthet(1,i,1,-1)=bthet(1,-i,1,1)
171 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
176 polthet(j,i)=polthet(j,-i)
179 gthet(j,i)=gthet(j,-i)
187 & 'Parameters of the virtual-bond valence angles:'
188 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
189 & ' ATHETA0 ',' A1 ',' A2 ',
192 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
193 & a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
195 write (iout,'(/a/9x,5a/79(1h-))')
196 & 'Parameters of the expression for sigma(theta_c):',
197 & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
198 & ' ALPH3 ',' SIGMA0C '
200 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
201 & (polthet(j,i),j=0,3),sigc0(i)
203 write (iout,'(/a/9x,5a/79(1h-))')
204 & 'Parameters of the second gaussian:',
205 & ' THETA0 ',' SIGMA0 ',' G1 ',
208 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
209 & sig0(i),(gthet(j,i),j=1,3)
213 & 'Parameters of the virtual-bond valence angles:'
214 write (iout,'(/a/9x,5a/79(1h-))')
215 & 'Coefficients of expansion',
216 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
217 & ' b1*10^1 ',' b2*10^1 '
219 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
220 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
221 & (10*bthet(j,i,1,1),j=1,2)
223 write (iout,'(/a/9x,5a/79(1h-))')
224 & 'Parameters of the expression for sigma(theta_c):',
225 & ' alpha0 ',' alph1 ',' alph2 ',
226 & ' alhp3 ',' sigma0c '
228 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
229 & (polthet(j,i),j=0,3),sigc0(i)
231 write (iout,'(/a/9x,5a/79(1h-))')
232 & 'Parameters of the second gaussian:',
233 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
236 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
237 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
242 IF (tor_mode.eq.0) THEN
244 C Read the parameters of Utheta determined from ab initio surfaces
245 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
247 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,
248 & ntheterm3,nsingle,ndouble
249 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
250 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
252 ithetyp(i)=-ithetyp(-i)
255 do i=-maxthetyp,maxthetyp
256 do j=-maxthetyp,maxthetyp
257 do k=-maxthetyp,maxthetyp
258 aa0thet(i,j,k,iblock)=0.0d0
260 aathet(l,i,j,k,iblock)=0.0d0
264 bbthet(m,l,i,j,k,iblock)=0.0d0
265 ccthet(m,l,i,j,k,iblock)=0.0d0
266 ddthet(m,l,i,j,k,iblock)=0.0d0
267 eethet(m,l,i,j,k,iblock)=0.0d0
273 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
274 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
282 c VAR:iblock means terminally blocking group 1=non-proline 2=proline
284 c VAR:ntethtyp is type of theta potentials type currently 0=glycine
285 c VAR:1=non-glicyne non-proline 2=proline
286 c VAR:negative values for D-aminoacid
288 do j=-nthetyp,nthetyp
289 do k=-nthetyp,nthetyp
290 read (ithep,'(6a)',end=111,err=111) res1
291 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
292 c VAR: aa0thet is variable describing the average value of Foureir
293 c VAR: expansion series
294 c VAR: aathet is foureir expansion in theta/2 angle for full formula
295 c VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
296 Condens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
297 read (ithep,*,end=111,err=111)
298 &(aathet(l,i,j,k,iblock),l=1,ntheterm)
299 read (ithep,*,end=111,err=111)
300 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
301 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
302 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
303 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),
305 read (ithep,*,end=111,err=111)
306 & (((ffthet(llll,lll,ll,i,j,k,iblock),
307 & ffthet(lll,llll,ll,i,j,k,iblock),
308 & ggthet(llll,lll,ll,i,j,k,iblock),
309 & ggthet(lll,llll,ll,i,j,k,iblock),
310 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
315 C For dummy ends assign glycine-type coefficients of theta-only terms; the
316 C coefficients of theta-and-gamma-dependent terms are zero.
317 C IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
318 C RECOMENTDED AFTER VERSION 3.3)
322 c aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
323 c aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
325 c aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
326 c aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
329 c aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
331 c aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
334 C AND COMMENT THE LOOPS BELOW
338 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
339 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
341 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
342 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
345 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
347 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
351 C Substitution for D aminoacids from symmetry.
354 do j=-nthetyp,nthetyp
355 do k=-nthetyp,nthetyp
356 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
358 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
362 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
363 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
364 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
365 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
371 ffthet(llll,lll,ll,i,j,k,iblock)=
372 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
373 ffthet(lll,llll,ll,i,j,k,iblock)=
374 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
375 ggthet(llll,lll,ll,i,j,k,iblock)=
376 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
377 ggthet(lll,llll,ll,i,j,k,iblock)=
378 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
387 C Control printout of the coefficients of virtual-bond-angle potentials
390 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
393 do j=-nthetyp,nthetyp
394 do k=-nthetyp,nthetyp
395 write (iout,'(//4a)')
396 & 'Type ',toronelet(i),toronelet(j),toronelet(k)
397 write (iout,'(//a,10x,a)') " l","a[l]"
398 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
399 write (iout,'(i2,1pe15.5)')
400 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
402 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))')
403 & "b",l,"c",l,"d",l,"e",l
405 write (iout,'(i2,4(1pe15.5))') m,
406 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
407 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
411 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
412 & "f+",l,"f-",l,"g+",l,"g-",l
415 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
416 & ffthet(n,m,l,i,j,k,iblock),
417 & ffthet(m,n,l,i,j,k,iblock),
418 & ggthet(n,m,l,i,j,k,iblock),
419 & ggthet(m,n,l,i,j,k,iblock)
432 C here will be the apropriate recalibrating for D-aminoacid
433 read (ithep,*,end=121,err=121) nthetyp
434 do i=-nthetyp+1,nthetyp-1
435 read (ithep,*,end=121,err=121) nbend_kcc_Tb(i)
436 do j=0,nbend_kcc_Tb(i)
437 read (ithep,*,end=121,err=121) ijunk,v1bend_chyb(j,i)
442 & "Parameters of the valence-only potentials"
443 do i=-nthetyp+1,nthetyp-1
444 write (iout,'(2a)') "Type ",toronelet(i)
445 do k=0,nbend_kcc_Tb(i)
446 write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
453 c write (2,*) "Start reading THETA_PDB",ithep_pdb
456 read (ithep_pdb,*,err=111,end=111)
457 & a0thet(i),(athet(j,i,1,1),j=1,2),
458 & (bthet(j,i,1,1),j=1,2)
459 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
460 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
461 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
465 athet(1,i,1,-1)=athet(1,i,1,1)
466 athet(2,i,1,-1)=athet(2,i,1,1)
467 bthet(1,i,1,-1)=-bthet(1,i,1,1)
468 bthet(2,i,1,-1)=-bthet(2,i,1,1)
469 athet(1,i,-1,1)=-athet(1,i,1,1)
470 athet(2,i,-1,1)=-athet(2,i,1,1)
471 bthet(1,i,-1,1)=bthet(1,i,1,1)
472 bthet(2,i,-1,1)=bthet(2,i,1,1)
476 athet(1,i,-1,-1)=athet(1,-i,1,1)
477 athet(2,i,-1,-1)=-athet(2,-i,1,1)
478 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
479 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
480 athet(1,i,-1,1)=athet(1,-i,1,1)
481 athet(2,i,-1,1)=-athet(2,-i,1,1)
482 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
483 bthet(2,i,-1,1)=bthet(2,-i,1,1)
484 athet(1,i,1,-1)=-athet(1,-i,1,1)
485 athet(2,i,1,-1)=athet(2,-i,1,1)
486 bthet(1,i,1,-1)=bthet(1,-i,1,1)
487 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
492 polthet(j,i)=polthet(j,-i)
495 gthet(j,i)=gthet(j,-i)
498 c write (2,*) "End reading THETA_PDB"
504 C Read the parameters of the probability distribution/energy expression
505 C of the side chains.
508 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
512 dsc_inv(i)=1.0D0/dsc(i)
523 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
524 & ((blower(k,l,1),l=1,k),k=1,3)
525 censc(1,1,-i)=censc(1,1,i)
526 censc(2,1,-i)=censc(2,1,i)
527 censc(3,1,-i)=-censc(3,1,i)
529 read (irotam,*,end=112,err=112) bsc(j,i)
530 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
531 & ((blower(k,l,j),l=1,k),k=1,3)
532 censc(1,j,-i)=censc(1,j,i)
533 censc(2,j,-i)=censc(2,j,i)
534 censc(3,j,-i)=-censc(3,j,i)
535 C BSC is amplitude of Gaussian
542 akl=akl+blower(k,m,j)*blower(l,m,j)
546 if (((k.eq.3).and.(l.ne.3))
547 & .or.((l.eq.3).and.(k.ne.3))) then
548 gaussc(k,l,j,-i)=-akl
549 gaussc(l,k,j,-i)=-akl
561 write (iout,'(/a)') 'Parameters of side-chain local geometry'
566 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
567 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
568 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
569 & 'log h',(bsc(j,i),j=1,nlobi)
570 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
571 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
573 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
574 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
577 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
578 write (iout,'(a,f10.4,4(16x,f10.4))')
579 & 'Center ',(bsc(j,i),j=1,nlobi)
580 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
589 C Read scrot parameters for potentials determined from all-atom AM1 calculations
590 C added by Urszula Kozlowska 07/11/2007
593 read (irotam,*,end=112,err=112)
595 read (irotam,*,end=112,err=112)
598 read(irotam,*,end=112,err=112) sc_parmin(j,i)
603 C Read the parameters of the probability distribution/energy expression
604 C of the side chains.
606 write (2,*) "Start reading ROTAM_PDB"
608 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
612 dsc_inv(i)=1.0D0/dsc(i)
623 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
624 & ((blower(k,l,1),l=1,k),k=1,3)
626 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
627 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
628 & ((blower(k,l,j),l=1,k),k=1,3)
635 akl=akl+blower(k,m,j)*blower(l,m,j)
645 write (2,*) "End reading ROTAM_PDB"
649 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
650 C interaction energy of the Gly, Ala, and Pro prototypes.
652 read (ifourier,*) nloctyp
653 SPLIT_FOURIERTOR = nloctyp.lt.0
654 nloctyp = iabs(nloctyp)
656 read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
657 read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
658 itype2loc(ntyp1)=nloctyp
659 iloctyp(nloctyp)=ntyp1
661 itype2loc(-i)=-itype2loc(i)
670 iloctyp(-i)=-iloctyp(i)
672 c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
673 c write (iout,*) "nloctyp",nloctyp,
674 c & " iloctyp",(iloctyp(i),i=0,nloctyp)
683 c write (iout,*) "NEWCORR",i
684 read (ifourier,*,end=115,err=115)
687 read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
690 c write (iout,*) "NEWCORR BNEW1"
691 c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
694 read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
697 c write (iout,*) "NEWCORR BNEW2"
698 c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
700 read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
701 read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
703 c write (iout,*) "NEWCORR CCNEW"
704 c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
706 read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
707 read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
709 c write (iout,*) "NEWCORR DDNEW"
710 c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
714 read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
718 c write (iout,*) "NEWCORR EENEW1"
719 c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
721 read (ifourier,*,end=115,err=115) e0new(ii,i)
723 c write (iout,*) (e0new(ii,i),ii=1,3)
725 c write (iout,*) "NEWCORR EENEW"
728 ccnew(ii,1,i)=ccnew(ii,1,i)/2
729 ccnew(ii,2,i)=ccnew(ii,2,i)/2
730 ddnew(ii,1,i)=ddnew(ii,1,i)/2
731 ddnew(ii,2,i)=ddnew(ii,2,i)/2
736 bnew1(ii,1,-i)= bnew1(ii,1,i)
737 bnew1(ii,2,-i)=-bnew1(ii,2,i)
738 bnew2(ii,1,-i)= bnew2(ii,1,i)
739 bnew2(ii,2,-i)=-bnew2(ii,2,i)
742 c ccnew(ii,1,i)=ccnew(ii,1,i)/2
743 c ccnew(ii,2,i)=ccnew(ii,2,i)/2
744 c ddnew(ii,1,i)=ddnew(ii,1,i)/2
745 c ddnew(ii,2,i)=ddnew(ii,2,i)/2
746 ccnew(ii,1,-i)=ccnew(ii,1,i)
747 ccnew(ii,2,-i)=-ccnew(ii,2,i)
748 ddnew(ii,1,-i)=ddnew(ii,1,i)
749 ddnew(ii,2,-i)=-ddnew(ii,2,i)
751 e0new(1,-i)= e0new(1,i)
752 e0new(2,-i)=-e0new(2,i)
753 e0new(3,-i)=-e0new(3,i)
755 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
756 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
757 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
758 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
762 write (iout,'(a)') "Coefficients of the multibody terms"
763 c do i=-nloctyp+1,nloctyp-1
764 do i=-nloctyp,nloctyp
765 write (iout,*) "Type: ",onelet(iloctyp(i))
766 write (iout,*) "Coefficients of the expansion of B1"
768 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
770 write (iout,*) "Coefficients of the expansion of B2"
772 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
774 write (iout,*) "Coefficients of the expansion of C"
775 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
776 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
777 write (iout,*) "Coefficients of the expansion of D"
778 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
779 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
780 write (iout,*) "Coefficients of the expansion of E"
781 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
784 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
789 IF (SPLIT_FOURIERTOR) THEN
791 c write (iout,*) "NEWCORR TOR",i
792 read (ifourier,*,end=115,err=115)
795 read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i)
798 c write (iout,*) "NEWCORR BNEW1 TOR"
799 c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
802 read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i)
805 c write (iout,*) "NEWCORR BNEW2 TOR"
806 c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
808 read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i)
809 read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i)
811 c write (iout,*) "NEWCORR CCNEW TOR"
812 c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
814 read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i)
815 read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i)
817 c write (iout,*) "NEWCORR DDNEW TOR"
818 c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
822 read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i)
826 c write (iout,*) "NEWCORR EENEW1 TOR"
827 c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
829 read (ifourier,*,end=115,err=115) e0newtor(ii,i)
831 c write (iout,*) (e0newtor(ii,i),ii=1,3)
833 c write (iout,*) "NEWCORR EENEW TOR"
836 ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
837 ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
838 ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
839 ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
844 bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
845 bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
846 bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
847 bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
850 ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
851 ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
852 ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
853 ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
855 e0newtor(1,-i)= e0newtor(1,i)
856 e0newtor(2,-i)=-e0newtor(2,i)
857 e0newtor(3,-i)=-e0newtor(3,i)
859 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
860 eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
861 eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
862 eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
867 & "Single-body coefficients of the torsional potentials"
868 do i=-nloctyp+1,nloctyp-1
869 write (iout,*) "Type: ",onelet(iloctyp(i))
870 write (iout,*) "Coefficients of the expansion of B1tor"
872 write (iout,'(3hB1(,i1,1h),3f10.5)')
873 & j,(bnew1tor(k,j,i),k=1,3)
875 write (iout,*) "Coefficients of the expansion of B2tor"
877 write (iout,'(3hB2(,i1,1h),3f10.5)')
878 & j,(bnew2tor(k,j,i),k=1,3)
880 write (iout,*) "Coefficients of the expansion of Ctor"
881 write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
882 write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
883 write (iout,*) "Coefficients of the expansion of Dtor"
884 write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
885 write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
886 write (iout,*) "Coefficients of the expansion of Etor"
887 write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
890 write (iout,'(1hE,2i1,2f10.5)')
891 & j,k,(eenewtor(l,j,k,i),l=1,2)
897 do i=-nloctyp+1,nloctyp-1
900 bnew1tor(ii,j,i)=bnew1(ii,j,i)
905 bnew2tor(ii,j,i)=bnew2(ii,j,i)
909 ccnewtor(ii,1,i)=ccnew(ii,1,i)
910 ccnewtor(ii,2,i)=ccnew(ii,2,i)
911 ddnewtor(ii,1,i)=ddnew(ii,1,i)
912 ddnewtor(ii,2,i)=ddnew(ii,2,i)
915 ENDIF !SPLIT_FOURIER_TOR
918 & write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
920 read (ifourier,*,end=115,err=115)
921 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
923 write (iout,*) 'Type ',onelet(iloctyp(i))
924 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
938 c B1tilde(1,i) = b(3)
939 c B1tilde(2,i) =-b(5)
940 c B1tilde(1,-i) =-b(3)
941 c B1tilde(2,-i) =b(5)
948 cc B1tilde(1,i) = b(3,i)
949 cc B1tilde(2,i) =-b(5,i)
950 C B1tilde(1,-i) =-b(3,i)
951 C B1tilde(2,-i) =b(5,i)
952 cc b1tilde(1,i)=0.0d0
953 cc b1tilde(2,i)=0.0d0
965 CCold(1,1,-i)= b(7,i)
966 CCold(2,2,-i)=-b(7,i)
967 CCold(2,1,-i)=-b(9,i)
968 CCold(1,2,-i)=-b(9,i)
973 c Ctilde(1,1,i)= CCold(1,1,i)
974 c Ctilde(1,2,i)= CCold(1,2,i)
975 c Ctilde(2,1,i)=-CCold(2,1,i)
976 c Ctilde(2,2,i)=-CCold(2,2,i)
978 c Ctilde(1,1,i)=0.0d0
979 c Ctilde(1,2,i)=0.0d0
980 c Ctilde(2,1,i)=0.0d0
981 c Ctilde(2,2,i)=0.0d0
986 DDold(1,1,-i)= b(6,i)
987 DDold(2,2,-i)=-b(6,i)
988 DDold(2,1,-i)=-b(8,i)
989 DDold(1,2,-i)=-b(8,i)
994 c Dtilde(1,1,i)= DD(1,1,i)
995 c Dtilde(1,2,i)= DD(1,2,i)
996 c Dtilde(2,1,i)=-DD(2,1,i)
997 c Dtilde(2,2,i)=-DD(2,2,i)
999 c Dtilde(1,1,i)=0.0d0
1000 c Dtilde(1,2,i)=0.0d0
1001 c Dtilde(2,1,i)=0.0d0
1002 c Dtilde(2,2,i)=0.0d0
1003 EEold(1,1,i)= b(10,i)+b(11,i)
1004 EEold(2,2,i)=-b(10,i)+b(11,i)
1005 EEold(2,1,i)= b(12,i)-b(13,i)
1006 EEold(1,2,i)= b(12,i)+b(13,i)
1007 EEold(1,1,-i)= b(10,i)+b(11,i)
1008 EEold(2,2,-i)=-b(10,i)+b(11,i)
1009 EEold(2,1,-i)=-b(12,i)+b(13,i)
1010 EEold(1,2,-i)=-b(12,i)-b(13,i)
1011 write(iout,*) "TU DOCHODZE"
1017 c ee(2,1,i)=ee(1,2,i)
1022 &"Coefficients of the cumulants (independent of valence angles)"
1023 do i=-nloctyp+1,nloctyp-1
1024 write (iout,*) 'Type ',onelet(iloctyp(i))
1026 write(iout,'(2f10.5)') B(3,i),B(5,i)
1028 write(iout,'(2f10.5)') B(2,i),B(4,i)
1031 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
1035 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
1039 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1047 C Read torsional parameters in old format
1049 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1050 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1051 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1056 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1062 write (iout,'(/a/)') 'Torsional constants:'
1065 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1066 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1072 C Read torsional parameters
1074 IF (TOR_MODE.eq.0) THEN
1076 read (itorp,*,end=113,err=113) ntortyp
1077 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1079 itype2loc(i)=itortyp(i)
1082 itype2loc(-i)=-itype2loc(i)
1084 itortyp(ntyp1)=ntortyp
1087 itortyp(i)=-itortyp(-i)
1089 write (iout,*) 'ntortyp',ntortyp
1091 do j=-ntortyp+1,ntortyp-1
1092 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
1094 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1095 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1098 do k=1,nterm(i,j,iblock)
1099 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
1101 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1102 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1103 v0ij=v0ij+si*v1(k,i,j,iblock)
1105 c write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
1106 c write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
1107 c &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
1109 do k=1,nlor(i,j,iblock)
1110 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
1111 & vlor2(k,i,j),vlor3(k,i,j)
1112 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1115 v0(-i,-j,iblock)=v0ij
1121 write (iout,'(/a/)') 'Parameters of the SCCOR potentials:'
1124 do j=-ntortyp+1,ntortyp-1
1125 write (iout,*) 'ityp',i,' jtyp',j
1126 write (iout,*) 'Fourier constants'
1127 do k=1,nterm(i,j,iblock)
1128 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
1131 write (iout,*) 'Lorenz constants'
1132 do k=1,nlor(i,j,iblock)
1133 write (iout,'(3(1pe15.5))')
1134 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1142 C 6/23/01 Read parameters for double torsionals
1146 do j=-ntortyp+1,ntortyp-1
1147 do k=-ntortyp+1,ntortyp-1
1148 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1149 c write (iout,*) "OK onelett",
1152 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
1153 & .or. t3.ne.toronelet(k)) then
1154 write (iout,*) "Error in double torsional parameter file",
1157 call MPI_Finalize(Ierror)
1159 stop "Error in double torsional parameter file"
1161 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
1162 & ntermd_2(i,j,k,iblock)
1163 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1164 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1165 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
1166 & ntermd_1(i,j,k,iblock))
1167 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
1168 & ntermd_1(i,j,k,iblock))
1169 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
1170 & ntermd_1(i,j,k,iblock))
1171 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
1172 & ntermd_1(i,j,k,iblock))
1173 C Martix of D parameters for one dimesional foureir series
1174 do l=1,ntermd_1(i,j,k,iblock)
1175 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1176 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1177 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1178 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1179 c write(iout,*) "whcodze" ,
1180 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1182 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
1183 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
1184 & v2s(m,l,i,j,k,iblock),
1185 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1186 C Martix of D parameters for two dimesional fourier series
1187 do l=1,ntermd_2(i,j,k,iblock)
1189 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1190 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1191 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1192 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1201 write (iout,*) 'Constants for double torsionals'
1204 do j=-ntortyp+1,ntortyp-1
1205 do k=-ntortyp+1,ntortyp-1
1206 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
1207 & ' nsingle',ntermd_1(i,j,k,iblock),
1208 & ' ndouble',ntermd_2(i,j,k,iblock)
1210 write (iout,*) 'Single angles:'
1211 do l=1,ntermd_1(i,j,k,iblock)
1212 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
1213 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
1214 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
1215 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1218 write (iout,*) 'Pairs of angles:'
1219 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1220 do l=1,ntermd_2(i,j,k,iblock)
1221 write (iout,'(i5,20f10.5)')
1222 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1225 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1226 do l=1,ntermd_2(i,j,k,iblock)
1227 write (iout,'(i5,20f10.5)')
1228 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
1229 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1238 ELSE IF (TOR_MODE.eq.1) THEN
1240 C read valence-torsional parameters
1241 read (itorp,*,end=121,err=121) ntortyp
1243 write (iout,*) "Valence-torsional parameters read in ntortyp",
1245 read (itorp,*,end=121,err=121) (itortyp(i),i=1,ntyp)
1246 write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
1249 itype2loc(i)=itortyp(i)
1253 itortyp(i)=-itortyp(-i)
1255 do i=-ntortyp+1,ntortyp-1
1256 do j=-ntortyp+1,ntortyp-1
1257 C first we read the cos and sin gamma parameters
1258 read (itorp,'(13x,a)',end=121,err=121) string
1259 write (iout,*) i,j,string
1260 read (itorp,*,end=121,err=121)
1261 & nterm_kcc(j,i),nterm_kcc_Tb(j,i)
1262 C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
1263 do k=1,nterm_kcc(j,i)
1264 do l=1,nterm_kcc_Tb(j,i)
1265 do ll=1,nterm_kcc_Tb(j,i)
1266 read (itorp,*,end=121,err=121) ii,jj,kk,
1267 & v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1275 c AL 4/8/16: Calculate coefficients from one-body parameters
1278 itortyp(i)=itype2loc(i)
1281 &"Val-tor parameters calculated from cumulant coefficients ntortyp"
1283 do i=-ntortyp+1,ntortyp-1
1284 do j=-ntortyp+1,ntortyp-1
1287 do k=1,nterm_kcc_Tb(j,i)
1288 do l=1,nterm_kcc_Tb(j,i)
1289 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)
1290 & +bnew1tor(k,2,i)*bnew2tor(l,2,j)
1291 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)
1292 & +bnew1tor(k,2,i)*bnew2tor(l,1,j)
1295 do k=1,nterm_kcc_Tb(j,i)
1296 do l=1,nterm_kcc_Tb(j,i)
1298 v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1299 & -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1300 v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1301 & +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1303 v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1304 & -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1305 v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1306 & +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1310 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)
1314 write (iout,*) "TOR_MODE>1 only with NEWCORR"
1319 if (tor_mode.gt.0 .and. lprint) then
1320 c Print valence-torsional parameters
1322 & "Parameters of the valence-torsional potentials"
1323 do i=-ntortyp+1,ntortyp-1
1324 do j=-ntortyp+1,ntortyp-1
1325 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
1326 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
1327 do k=1,nterm_kcc(j,i)
1328 do l=1,nterm_kcc_Tb(j,i)
1329 do ll=1,nterm_kcc_Tb(j,i)
1330 write (iout,'(3i5,2f15.4)')
1331 & k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1340 C Read of Side-chain backbone correlation parameters
1341 C Modified 11 May 2012 by Adasko
1344 read (isccor,*,end=119,err=119) nsccortyp
1346 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1348 isccortyp(i)=-isccortyp(-i)
1350 iscprol=isccortyp(20)
1351 c write (iout,*) 'ntortyp',ntortyp
1353 cc maxinter is maximum interaction sites
1357 read (isccor,*,end=119,err=119)
1358 &nterm_sccor(i,j),nlor_sccor(i,j)
1364 nterm_sccor(-i,j)=nterm_sccor(i,j)
1365 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1366 nterm_sccor(i,-j)=nterm_sccor(i,j)
1367 do k=1,nterm_sccor(i,j)
1368 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
1370 if (j.eq.iscprol) then
1371 if (i.eq.isccortyp(10)) then
1372 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1373 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1375 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
1376 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1377 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
1378 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1379 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1380 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1381 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1382 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1385 if (i.eq.isccortyp(10)) then
1386 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1387 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1389 if (j.eq.isccortyp(10)) then
1390 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1391 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1393 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1394 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1395 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1396 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1397 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1398 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1402 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1403 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1404 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1405 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1408 do k=1,nlor_sccor(i,j)
1409 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
1410 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1411 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
1412 &(1+vlor3sccor(k,i,j)**2)
1414 v0sccor(l,i,j)=v0ijsccor
1415 v0sccor(l,-i,j)=v0ijsccor1
1416 v0sccor(l,i,-j)=v0ijsccor2
1417 v0sccor(l,-i,-j)=v0ijsccor3
1423 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1424 c write (iout,*) 'ntortyp',ntortyp
1426 cc maxinter is maximum interaction sites
1430 read (isccor,*,end=119,err=119)
1431 & nterm_sccor(i,j),nlor_sccor(i,j)
1435 do k=1,nterm_sccor(i,j)
1436 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
1438 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1441 do k=1,nlor_sccor(i,j)
1442 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
1443 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1444 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
1445 &(1+vlor3sccor(k,i,j)**2)
1447 v0sccor(l,i,j)=v0ijsccor
1455 write (iout,'(/a/)') 'Torsional constants:'
1459 write (iout,*) 'ityp',i,' jtyp',j
1460 write (iout,*) 'Fourier constants'
1461 do k=1,nterm_sccor(i,j)
1462 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1464 write (iout,*) 'Lorenz constants'
1465 do k=1,nlor_sccor(i,j)
1466 write (iout,'(3(1pe15.5))')
1467 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1475 C Read electrostatic-interaction parameters
1479 write (iout,'(/a)') 'Electrostatic interaction constants:'
1480 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1481 & 'IT','JT','APP','BPP','AEL6','AEL3'
1483 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1484 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1485 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1486 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1491 app (i,j)=epp(i,j)*rri*rri
1492 bpp (i,j)=-2.0D0*epp(i,j)*rri
1493 ael6(i,j)=elpp6(i,j)*4.2D0**6
1494 ael3(i,j)=elpp3(i,j)*4.2D0**3
1496 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1497 & ael6(i,j),ael3(i,j)
1502 C Read side-chain interaction parameters.
1504 read (isidep,*,end=117,err=117) ipot,expon
1505 if (ipot.lt.1 .or. ipot.gt.5) then
1506 write (iout,'(2a)') 'Error while reading SC interaction',
1507 & 'potential file - unknown potential type.'
1509 call MPI_Finalize(Ierror)
1515 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1516 & ', exponents are ',expon,2*expon
1517 goto (10,20,30,30,40) ipot
1518 C----------------------- LJ potential ---------------------------------
1519 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1520 & (sigma0(i),i=1,ntyp)
1522 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1523 write (iout,'(a/)') 'The epsilon array:'
1524 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1525 write (iout,'(/a)') 'One-body parameters:'
1526 write (iout,'(a,4x,a)') 'residue','sigma'
1527 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1530 C----------------------- LJK potential --------------------------------
1531 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1532 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1534 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1535 write (iout,'(a/)') 'The epsilon array:'
1536 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1537 write (iout,'(/a)') 'One-body parameters:'
1538 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1539 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1543 C---------------------- GB or BP potential -----------------------------
1545 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1547 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
1548 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
1549 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
1550 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
1551 C now we start reading lipid
1553 read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
1554 if (lprint) write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp)
1556 C print *,"WARNING!!"
1558 C epslip(i,j)=epslip(i,j)+0.05d0
1561 if (lprint) write(iout,*) epslip(1,1),"OK?"
1562 C For the GB potential convert sigma'**2 into chi'
1565 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1569 write (iout,'(/a/)') 'Parameters of the BP potential:'
1570 write (iout,'(a/)') 'The epsilon array:'
1571 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1572 write (iout,'(/a)') 'One-body parameters:'
1573 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1575 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1576 & chip(i),alp(i),i=1,ntyp)
1579 C--------------------- GBV potential -----------------------------------
1580 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1581 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1582 & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1584 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1585 write (iout,'(a/)') 'The epsilon array:'
1586 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1587 write (iout,'(/a)') 'One-body parameters:'
1588 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1589 & 's||/s_|_^2',' chip ',' alph '
1590 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1591 & sigii(i),chip(i),alp(i),i=1,ntyp)
1595 C-----------------------------------------------------------------------
1596 C Calculate the "working" parameters of SC interactions.
1600 epslip(i,j)=epslip(j,i)
1605 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1606 sigma(j,i)=sigma(i,j)
1607 rs0(i,j)=dwa16*sigma(i,j)
1611 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1612 & 'Working parameters of the SC interactions:',
1613 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1618 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1627 sigeps=dsign(1.0D0,epsij)
1629 aa_aq(i,j)=epsij*rrij*rrij
1630 bb_aq(i,j)=-sigeps*epsij*rrij
1631 aa_aq(j,i)=aa_aq(i,j)
1632 bb_aq(j,i)=bb_aq(i,j)
1633 epsijlip=epslip(i,j)
1634 sigeps=dsign(1.0D0,epsijlip)
1635 epsijlip=dabs(epsijlip)
1636 aa_lip(i,j)=epsijlip*rrij*rrij
1637 bb_lip(i,j)=-sigeps*epsijlip*rrij
1638 aa_lip(j,i)=aa_lip(i,j)
1639 bb_lip(j,i)=bb_lip(i,j)
1641 sigt1sq=sigma0(i)**2
1642 sigt2sq=sigma0(j)**2
1645 ratsig1=sigt2sq/sigt1sq
1646 ratsig2=1.0D0/ratsig1
1647 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1648 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1649 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1653 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1654 sigmaii(i,j)=rsum_max
1655 sigmaii(j,i)=rsum_max
1657 c sigmaii(i,j)=r0(i,j)
1658 c sigmaii(j,i)=r0(i,j)
1660 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1661 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1662 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1663 augm(i,j)=epsij*r_augm**(2*expon)
1664 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1671 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1672 & restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1673 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1678 write(iout,*) "tube param"
1679 read(itube,*) epspeptube,sigmapeptube
1680 sigmapeptube=sigmapeptube**6
1681 sigeps=dsign(1.0D0,epspeptube)
1682 epspeptube=dabs(epspeptube)
1683 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
1684 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
1685 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
1687 read(itube,*) epssctube,sigmasctube
1688 sigmasctube=sigmasctube**6
1689 sigeps=dsign(1.0D0,epssctube)
1690 epssctube=dabs(epssctube)
1691 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
1692 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
1693 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
1699 C Define the SC-p interaction constants (hard-coded; old style)
1702 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1704 c aad(i,1)=0.3D0*4.0D0**12
1705 C Following line for constants currently implemented
1706 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1707 aad(i,1)=1.5D0*4.0D0**12
1708 c aad(i,1)=0.17D0*5.6D0**12
1710 C "Soft" SC-p repulsion
1712 C Following line for constants currently implemented
1713 c aad(i,1)=0.3D0*4.0D0**6
1714 C "Hard" SC-p repulsion
1715 bad(i,1)=3.0D0*4.0D0**6
1716 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1725 C 8/9/01 Read the SC-p interaction constants from file
1728 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1731 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1732 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1733 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1734 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1738 write (iout,*) "Parameters of SC-p interactions:"
1740 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1741 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1747 C Define the constants of the disulfide bridge
1751 c Old arbitrary potential - commented out.
1756 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1757 c energy surface of diethyl disulfide.
1758 c A. Liwo and U. Kozlowska, 11/24/03
1774 C if(me.eq.king) then
1775 C write (iout,'(/a)') "Disulfide bridge parameters:"
1776 C write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1777 C write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1778 C write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1779 C write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1782 if (shield_mode.gt.0) then
1783 C VSolvSphere the volume of solving sphere
1784 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1785 C there will be no distinction between proline peptide group and normal peptide
1786 C group in case of shielding parameters
1787 c write (iout,*) "rpp(1,1)",rpp(1,1)," pi",pi
1788 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1789 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1790 write (iout,*) "VSolvSphere",VSolvSphere,"VSolvSphere_div",
1792 C long axis of side chain
1794 long_r_sidechain(i)=vbldsc0(1,i)
1795 short_r_sidechain(i)=sigma0(i)
1800 111 write (iout,*) "Error reading bending energy parameters."
1802 112 write (iout,*) "Error reading rotamer energy parameters."
1804 113 write (iout,*) "Error reading torsional energy parameters."
1806 114 write (iout,*) "Error reading double torsional energy parameters."
1809 & "Error reading cumulant (multibody energy) parameters."
1811 116 write (iout,*) "Error reading electrostatic energy parameters."
1813 1161 write (iout,*) "Error reading electrostatic energy parameters.Lip"
1815 117 write (iout,*) "Error reading side chain interaction parameters."
1817 118 write (iout,*) "Error reading SCp interaction parameters."
1819 119 write (iout,*) "Error reading SCCOR parameters"
1821 121 write (iout,*) "Error in Czybyshev parameters"
1824 call MPI_Finalize(Ierror)
1831 subroutine getenv_loc(var, val)
1832 character(*) var, val
1835 character(2000) line
1838 open (196,file='env',status='old',readonly,shared)
1840 c write(*,*)'looking for ',var
1841 10 read(196,*,err=11,end=11)line
1842 iread=index(line,var)
1843 c write(*,*)iread,' ',var,' ',line
1844 if (iread.eq.0) go to 10
1845 c write(*,*)'---> ',line
1851 iread=iread+ilen(var)+1
1852 read (line(iread:),*,err=12,end=12) val
1853 c write(*,*)'OK: ',var,' = ',val
1859 #elif (defined CRAY)
1860 integer lennam,lenval,ierror
1862 c getenv using a POSIX call, useful on the T3D
1863 c Sept 1996, comment out error check on advice of H. Pritchard
1866 if(lennam.le.0) stop '--error calling getenv--'
1867 call pxfgetenv(var,lennam,val,lenval,ierror)
1868 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1870 call getenv(var,val)
1872 C set the variables used for shielding effect
1873 C if (shield_mode.gt.0) then
1874 C VSolvSphere the volume of solving sphere
1876 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1877 C there will be no distinction between proline peptide group and normal peptide
1878 C group in case of shielding parameters
1879 C VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1880 C VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1881 C long axis of side chain
1883 C long_r_sidechain(i)=vbldsc0(1,i)
1884 C short_r_sidechain(i)=sigma0(i)
1886 C lets set the buffor value
1891 c---------------------------------------------------------------------------
1892 subroutine weightread
1894 C Read the parameters of the probability distributions of the virtual-bond
1895 C valence angles and the side chains and energy parameters.
1897 C Important! Energy-term weights ARE NOT read here; they are read from the
1898 C main input file instead, because NO defaults have yet been set for these
1902 include 'DIMENSIONS'
1907 include 'COMMON.IOUNITS'
1908 include 'COMMON.CHAIN'
1909 include 'COMMON.INTERACT'
1910 include 'COMMON.GEO'
1911 include 'COMMON.LOCAL'
1912 include 'COMMON.TORSION'
1913 include 'COMMON.SCCOR'
1914 include 'COMMON.SCROT'
1915 include 'COMMON.FFIELD'
1916 include 'COMMON.NAMES'
1917 include 'COMMON.SBRIDGE'
1919 include 'COMMON.SETUP'
1920 include 'COMMON.CONTROL'
1921 include 'COMMON.SHIELD'
1922 character*1000 weightcard
1924 double precision scalscp,wlong
1926 c READ energy-term weights
1928 call card_concat(weightcard)
1929 call reada(weightcard,'WLONG',wlong,1.0D0)
1930 call reada(weightcard,'WSC',wsc,wlong)
1931 call reada(weightcard,'WSCP',wscp,wlong)
1932 call reada(weightcard,'WELEC',welec,1.0D0)
1933 call reada(weightcard,'WVDWPP',wvdwpp,welec)
1934 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
1935 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
1936 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
1937 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
1938 call reada(weightcard,'WTURN3',wturn3,1.0D0)
1939 call reada(weightcard,'WTURN4',wturn4,1.0D0)
1940 call reada(weightcard,'WTURN6',wturn6,1.0D0)
1941 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
1942 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
1943 call reada(weightcard,'WBOND',wbond,1.0D0)
1944 call reada(weightcard,'WTOR',wtor,1.0D0)
1945 call reada(weightcard,'WTORD',wtor_d,1.0D0)
1946 call reada(weightcard,'WANG',wang,1.0D0)
1947 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
1948 call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
1949 call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
1950 call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
1951 call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
1952 call reada(weightcard,'SCAL14',scal14,0.4D0)
1953 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
1954 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
1955 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
1956 call reada(weightcard,'TEMP0',temp0,300.0d0)
1957 call reada(weightcard,'WSHIELD',wshield,1.0d0)
1958 call reada(weightcard,'WLT',wliptran,0.0D0)
1959 call reada(weightcard,'WTUBE',wtube,1.0D0)
1960 call reada(weightcard,'WSAXS',wsaxs,1.0D0)
1961 if (index(weightcard,'SOFT').gt.0) ipot=6
1962 C 12/1/95 Added weight for the multi-body term WCORR
1963 call reada(weightcard,'WCORRH',wcorr,1.0D0)
1964 if (wcorr4.gt.0.0d0) wcorr=wcorr4
1985 weights(28)=wdfa_dist
1986 weights(29)=wdfa_tor
1987 weights(30)=wdfa_nei
1988 weights(31)=wdfa_beta
1990 aad(i,1)=scalscp*aad(i,1)
1991 aad(i,2)=scalscp*aad(i,2)
1992 bad(i,1)=scalscp*bad(i,1)
1993 bad(i,2)=scalscp*bad(i,2)
1995 if(me.eq.king.or..not.out1file)
1996 & write (iout,100) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
1997 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
1999 100 format (/'Energy-term weights (unscaled):'//
2000 & 'WSCC= ',f10.6,' (SC-SC)'/
2001 & 'WSCP= ',f10.6,' (SC-p)'/
2002 & 'WELEC= ',f10.6,' (p-p electr)'/
2003 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
2004 & 'WBOND= ',f10.6,' (stretching)'/
2005 & 'WANG= ',f10.6,' (bending)'/
2006 & 'WSCLOC= ',f10.6,' (SC local)'/
2007 & 'WTOR= ',f10.6,' (torsional)'/
2008 & 'WTORD= ',f10.6,' (double torsional)'/
2009 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
2010 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
2011 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
2012 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
2013 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
2014 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
2015 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
2016 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
2017 & 'WTURN6= ',f10.6,' (turns, 6th order)')
2018 if(me.eq.king.or..not.out1file)then
2019 if (wcorr4.gt.0.0d0) then
2020 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
2021 & 'between contact pairs of peptide groups'
2022 write (iout,'(2(a,f5.3/))')
2023 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
2024 & 'Range of quenching the correlation terms:',2*delt_corr
2025 else if (wcorr.gt.0.0d0) then
2026 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
2027 & 'between contact pairs of peptide groups'
2029 write (iout,'(a,f8.3)')
2030 & 'Scaling factor of 1,4 SC-p interactions:',scal14
2031 write (iout,'(a,f8.3)')
2032 & 'General scaling factor of SC-p interactions:',scalscp
2034 r0_corr=cutoff_corr-delt_corr
2036 call rescale_weights(t_bath)
2037 if(me.eq.king.or..not.out1file)
2038 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
2039 & wtor_d,wstrain,wel_loc,
2041 & wcorr,wcorr5,wcorr6,
2048 22 format (/'Energy-term weights (scaled):'//
2049 & 'WSCC= ',f10.6,' (SC-SC)'/
2050 & 'WSCP= ',f10.6,' (SC-p)'/
2051 & 'WELEC= ',f10.6,' (p-p electr)'/
2052 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
2053 & 'WBOND= ',f10.6,' (stretching)'/
2054 & 'WANG= ',f10.6,' (bending)'/
2055 & 'WSCLOC= ',f10.6,' (SC local)'/
2056 & 'WTOR= ',f10.6,' (torsional)'/
2057 & 'WTORD= ',f10.6,' (double torsional)'/
2058 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
2059 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
2061 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
2062 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
2063 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
2065 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
2066 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
2067 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
2069 & 'WTURN6= ',f10.6,' (turns, 6th order)'
2072 if(me.eq.king.or..not.out1file)
2073 & write (iout,*) "Reference temperature for weights calculation:",
2076 write (iout,'(/a/)') "DFA pseudopotential parameters:"
2077 write (iout,'(a,f10.6,a)')
2078 & "WDFAD= ",wdfa_dist," (distance)",
2079 & "WDFAT= ",wdfa_tor," (backbone angles)",
2080 & "WDFAN= ",wdfa_nei," (neighbors)",
2081 & "WDFAB= ",wdfa_beta," (beta structure)"
2083 call reada(weightcard,"D0CM",d0cm,3.78d0)
2084 call reada(weightcard,"AKCM",akcm,15.1d0)
2085 call reada(weightcard,"AKTH",akth,11.0d0)
2086 call reada(weightcard,"AKCT",akct,12.0d0)
2087 call reada(weightcard,"V1SS",v1ss,-1.08d0)
2088 call reada(weightcard,"V2SS",v2ss,7.61d0)
2089 call reada(weightcard,"V3SS",v3ss,13.7d0)
2090 call reada(weightcard,"EBR",ebr,-5.50D0)
2091 call reada(weightcard,"ATRISS",atriss,0.301D0)
2092 call reada(weightcard,"BTRISS",btriss,0.021D0)
2093 call reada(weightcard,"CTRISS",ctriss,1.001D0)
2094 call reada(weightcard,"DTRISS",dtriss,1.001D0)
2095 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
2097 dyn_ss_mask(i)=.false.
2101 dyn_ssbond_ij(i,j)=1.0d300
2104 call reada(weightcard,"HT",Ht,0.0D0)
2106 ss_depth=ebr/wsc-0.25*eps(1,1)
2107 Ht=Ht/wsc-0.25*eps(1,1)
2108 akcm=akcm*wstrain/wsc
2109 akth=akth*wstrain/wsc
2110 akct=akct*wstrain/wsc
2111 v1ss=v1ss*wstrain/wsc
2112 v2ss=v2ss*wstrain/wsc
2113 v3ss=v3ss*wstrain/wsc
2115 if (wstrain.ne.0.0) then
2116 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
2121 write (iout,*) "wshield,", wshield
2122 if(me.eq.king.or..not.out1file) then
2123 write (iout,*) "Parameters of the SS-bond potential:"
2124 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
2126 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
2127 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
2128 write (iout,*)" HT",Ht
2129 write (iout,*) "Parameters of the 'trisulfide' potential"
2130 write (iout,*) "ATRISS=", atriss
2131 write (iout,*) "BTRISS=", btriss
2132 write (iout,*) "CTRISS=", ctriss
2133 write (iout,*) "DTRISS=", dtriss
2134 print *,'indpdb=',indpdb,' pdbref=',pdbref