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
10 implicit real*8 (a-h,o-z)
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'
28 include 'COMMON.SETUP'
29 include 'COMMON.CONTROL'
30 include 'COMMON.SHIELD'
32 character*1 onelett(4) /"G","A","P","D"/
33 character*1 toronelet(-2:2) /"p","a","G","A","P"/
35 dimension blower(3,3,maxlob)
38 character*3 lancuch,ucase
40 C For printing parameters after they are read set the following in the UNRES
43 C setenv PRINT_PARM YES
45 C To print parameters in LaTeX format rather than as ASCII tables:
49 call getenv_loc("PRINT_PARM",lancuch)
50 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
51 call getenv_loc("LATEX",lancuch)
52 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
54 dwa16=2.0d0**(1.0d0/6.0d0)
56 C Assign virtual-bond length
61 c Read the virtual-bond parameters, masses, and moments of inertia
62 c and Stokes radii of the peptide group and side chains
65 read (ibond,*) vbldp0,vbldpdum,akp,mp,ip,pstok
68 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
73 dsc_inv(i)=1.0D0/dsc(i)
77 read (ibond,*) junk,vbldp0,vbldpdum,akp,rjunk,mp,ip,pstok
79 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
80 & j=1,nbondterm(i)),msc(i),isc(i),restok(i)
85 dsc_inv(i)=1.0D0/dsc(i)
90 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
91 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',
93 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
95 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
96 & vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
98 write (iout,'(13x,3f10.5)')
99 & vbldsc0(j,i),aksc(j,i),abond0(j,i)
103 C reading lipid parameters
105 write (iout,*) "iliptranpar",iliptranpar
108 read(iliptranpar,*) pepliptran
110 read(iliptranpar,*) liptranene(i)
115 C Read the parameters of the probability distribution/energy expression
116 C of the virtual-bond valence angles theta
119 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),
120 & (bthet(j,i,1,1),j=1,2)
121 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
122 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
123 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
127 athet(1,i,1,-1)=athet(1,i,1,1)
128 athet(2,i,1,-1)=athet(2,i,1,1)
129 bthet(1,i,1,-1)=-bthet(1,i,1,1)
130 bthet(2,i,1,-1)=-bthet(2,i,1,1)
131 athet(1,i,-1,1)=-athet(1,i,1,1)
132 athet(2,i,-1,1)=-athet(2,i,1,1)
133 bthet(1,i,-1,1)=bthet(1,i,1,1)
134 bthet(2,i,-1,1)=bthet(2,i,1,1)
138 athet(1,i,-1,-1)=athet(1,-i,1,1)
139 athet(2,i,-1,-1)=-athet(2,-i,1,1)
140 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
141 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
142 athet(1,i,-1,1)=athet(1,-i,1,1)
143 athet(2,i,-1,1)=-athet(2,-i,1,1)
144 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
145 bthet(2,i,-1,1)=bthet(2,-i,1,1)
146 athet(1,i,1,-1)=-athet(1,-i,1,1)
147 athet(2,i,1,-1)=athet(2,-i,1,1)
148 bthet(1,i,1,-1)=bthet(1,-i,1,1)
149 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
154 polthet(j,i)=polthet(j,-i)
157 gthet(j,i)=gthet(j,-i)
165 & 'Parameters of the virtual-bond valence angles:'
166 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
167 & ' ATHETA0 ',' A1 ',' A2 ',
170 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
171 & a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
173 write (iout,'(/a/9x,5a/79(1h-))')
174 & 'Parameters of the expression for sigma(theta_c):',
175 & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
176 & ' ALPH3 ',' SIGMA0C '
178 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
179 & (polthet(j,i),j=0,3),sigc0(i)
181 write (iout,'(/a/9x,5a/79(1h-))')
182 & 'Parameters of the second gaussian:',
183 & ' THETA0 ',' SIGMA0 ',' G1 ',
186 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
187 & sig0(i),(gthet(j,i),j=1,3)
191 & 'Parameters of the virtual-bond valence angles:'
192 write (iout,'(/a/9x,5a/79(1h-))')
193 & 'Coefficients of expansion',
194 & ' theta0 ',' a1*10^2 ',' a2*10^2 ',
195 & ' b1*10^1 ',' b2*10^1 '
197 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
198 & a0thet(i),(100*athet(j,i,1,1),j=1,2),
199 & (10*bthet(j,i,1,1),j=1,2)
201 write (iout,'(/a/9x,5a/79(1h-))')
202 & 'Parameters of the expression for sigma(theta_c):',
203 & ' alpha0 ',' alph1 ',' alph2 ',
204 & ' alhp3 ',' sigma0c '
206 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
207 & (polthet(j,i),j=0,3),sigc0(i)
209 write (iout,'(/a/9x,5a/79(1h-))')
210 & 'Parameters of the second gaussian:',
211 & ' theta0 ',' sigma0*10^2 ',' G1*10^-1',
214 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
215 & 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
220 IF (tor_mode.eq.0) THEN
222 C Read the parameters of Utheta determined from ab initio surfaces
223 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
225 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,
226 & ntheterm3,nsingle,ndouble
227 write (iout,*) "ithep",ithep
229 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
230 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
232 ithetyp(i)=-ithetyp(-i)
235 do i=-maxthetyp,maxthetyp
236 do j=-maxthetyp,maxthetyp
237 do k=-maxthetyp,maxthetyp
238 aa0thet(i,j,k,iblock)=0.0d0
240 aathet(l,i,j,k,iblock)=0.0d0
244 bbthet(m,l,i,j,k,iblock)=0.0d0
245 ccthet(m,l,i,j,k,iblock)=0.0d0
246 ddthet(m,l,i,j,k,iblock)=0.0d0
247 eethet(m,l,i,j,k,iblock)=0.0d0
253 ffthet(mm,m,l,i,j,k,iblock)=0.0d0
254 ggthet(mm,m,l,i,j,k,iblock)=0.0d0
262 c VAR:iblock means terminally blocking group 1=non-proline 2=proline
264 c VAR:ntethtyp is type of theta potentials type currently 0=glycine
265 c VAR:1=non-glicyne non-proline 2=proline
266 c VAR:negative values for D-aminoacid
268 do j=-nthetyp,nthetyp
269 do k=-nthetyp,nthetyp
270 read (ithep,'(6a)',end=111,err=111) res1
271 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
272 c VAR: aa0thet is variable describing the average value of Foureir
273 c VAR: expansion series
274 c VAR: aathet is foureir expansion in theta/2 angle for full formula
275 c VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
276 Condens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
277 read (ithep,*,end=111,err=111)
278 &(aathet(l,i,j,k,iblock),l=1,ntheterm)
279 read (ithep,*,end=111,err=111)
280 & ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
281 & (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
282 & (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
283 & (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),
285 read (ithep,*,end=111,err=111)
286 & (((ffthet(llll,lll,ll,i,j,k,iblock),
287 & ffthet(lll,llll,ll,i,j,k,iblock),
288 & ggthet(llll,lll,ll,i,j,k,iblock),
289 & ggthet(lll,llll,ll,i,j,k,iblock),
290 & llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
295 C For dummy ends assign glycine-type coefficients of theta-only terms; the
296 C coefficients of theta-and-gamma-dependent terms are zero.
297 C IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
298 C RECOMENTDED AFTER VERSION 3.3)
302 c aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
303 c aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
305 c aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
306 c aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
309 c aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
311 c aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
314 C AND COMMENT THE LOOPS BELOW
318 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
319 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
321 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
322 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
325 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
327 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
331 C Substitution for D aminoacids from symmetry.
334 do j=-nthetyp,nthetyp
335 do k=-nthetyp,nthetyp
336 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
338 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
342 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
343 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
344 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
345 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
351 ffthet(llll,lll,ll,i,j,k,iblock)=
352 & ffthet(llll,lll,ll,-i,-j,-k,iblock)
353 ffthet(lll,llll,ll,i,j,k,iblock)=
354 & ffthet(lll,llll,ll,-i,-j,-k,iblock)
355 ggthet(llll,lll,ll,i,j,k,iblock)=
356 & -ggthet(llll,lll,ll,-i,-j,-k,iblock)
357 ggthet(lll,llll,ll,i,j,k,iblock)=
358 & -ggthet(lll,llll,ll,-i,-j,-k,iblock)
367 C Control printout of the coefficients of virtual-bond-angle potentials
370 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
373 do j=-nthetyp,nthetyp
374 do k=-nthetyp,nthetyp
375 write (iout,'(//4a)')
376 & 'Type ',toronelet(i),toronelet(j),toronelet(k)
377 write (iout,'(//a,10x,a)') " l","a[l]"
378 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
379 write (iout,'(i2,1pe15.5)')
380 & (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
382 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))')
383 & "b",l,"c",l,"d",l,"e",l
385 write (iout,'(i2,4(1pe15.5))') m,
386 & bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
387 & ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
391 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
392 & "f+",l,"f-",l,"g+",l,"g-",l
395 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
396 & ffthet(n,m,l,i,j,k,iblock),
397 & ffthet(m,n,l,i,j,k,iblock),
398 & ggthet(n,m,l,i,j,k,iblock),
399 & ggthet(m,n,l,i,j,k,iblock)
412 C here will be the apropriate recalibrating for D-aminoacid
413 read (ithep,*,end=121,err=121) nthetyp
414 do i=-nthetyp+1,nthetyp-1
415 read (ithep,*,end=121,err=121) nbend_kcc_Tb(i)
416 do j=0,nbend_kcc_Tb(i)
417 read (ithep,*,end=121,err=121) ijunk,v1bend_chyb(j,i)
422 & "Parameters of the valence-only potentials"
423 do i=-nthetyp+1,nthetyp-1
424 write (iout,'(2a)') "Type ",toronelet(i)
425 do k=0,nbend_kcc_Tb(i)
426 write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
433 c write (2,*) "Start reading THETA_PDB",ithep_pdb
436 read (ithep_pdb,*,err=111,end=111)
437 & a0thet(i),(athet(j,i,1,1),j=1,2),
438 & (bthet(j,i,1,1),j=1,2)
439 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
440 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
441 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
445 athet(1,i,1,-1)=athet(1,i,1,1)
446 athet(2,i,1,-1)=athet(2,i,1,1)
447 bthet(1,i,1,-1)=-bthet(1,i,1,1)
448 bthet(2,i,1,-1)=-bthet(2,i,1,1)
449 athet(1,i,-1,1)=-athet(1,i,1,1)
450 athet(2,i,-1,1)=-athet(2,i,1,1)
451 bthet(1,i,-1,1)=bthet(1,i,1,1)
452 bthet(2,i,-1,1)=bthet(2,i,1,1)
456 athet(1,i,-1,-1)=athet(1,-i,1,1)
457 athet(2,i,-1,-1)=-athet(2,-i,1,1)
458 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
459 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
460 athet(1,i,-1,1)=athet(1,-i,1,1)
461 athet(2,i,-1,1)=-athet(2,-i,1,1)
462 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
463 bthet(2,i,-1,1)=bthet(2,-i,1,1)
464 athet(1,i,1,-1)=-athet(1,-i,1,1)
465 athet(2,i,1,-1)=athet(2,-i,1,1)
466 bthet(1,i,1,-1)=bthet(1,-i,1,1)
467 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
472 polthet(j,i)=polthet(j,-i)
475 gthet(j,i)=gthet(j,-i)
478 c write (2,*) "End reading THETA_PDB"
484 C Read the parameters of the probability distribution/energy expression
485 C of the side chains.
488 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
492 dsc_inv(i)=1.0D0/dsc(i)
503 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
504 & ((blower(k,l,1),l=1,k),k=1,3)
505 censc(1,1,-i)=censc(1,1,i)
506 censc(2,1,-i)=censc(2,1,i)
507 censc(3,1,-i)=-censc(3,1,i)
509 read (irotam,*,end=112,err=112) bsc(j,i)
510 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
511 & ((blower(k,l,j),l=1,k),k=1,3)
512 censc(1,j,-i)=censc(1,j,i)
513 censc(2,j,-i)=censc(2,j,i)
514 censc(3,j,-i)=-censc(3,j,i)
515 C BSC is amplitude of Gaussian
522 akl=akl+blower(k,m,j)*blower(l,m,j)
526 if (((k.eq.3).and.(l.ne.3))
527 & .or.((l.eq.3).and.(k.ne.3))) then
528 gaussc(k,l,j,-i)=-akl
529 gaussc(l,k,j,-i)=-akl
541 write (iout,'(/a)') 'Parameters of side-chain local geometry'
546 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
547 & ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
548 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
549 & 'log h',(bsc(j,i),j=1,nlobi)
550 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
551 & 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
553 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
554 & ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
557 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
558 write (iout,'(a,f10.4,4(16x,f10.4))')
559 & 'Center ',(bsc(j,i),j=1,nlobi)
560 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
569 C Read scrot parameters for potentials determined from all-atom AM1 calculations
570 C added by Urszula Kozlowska 07/11/2007
573 read (irotam,*,end=112,err=112)
575 read (irotam,*,end=112,err=112)
578 read(irotam,*,end=112,err=112) sc_parmin(j,i)
583 C Read the parameters of the probability distribution/energy expression
584 C of the side chains.
586 write (2,*) "Start reading ROTAM_PDB"
588 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
592 dsc_inv(i)=1.0D0/dsc(i)
603 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
604 & ((blower(k,l,1),l=1,k),k=1,3)
606 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
607 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
608 & ((blower(k,l,j),l=1,k),k=1,3)
615 akl=akl+blower(k,m,j)*blower(l,m,j)
625 write (2,*) "End reading ROTAM_PDB"
629 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
630 C interaction energy of the Gly, Ala, and Pro prototypes.
632 read (ifourier,*) nloctyp
633 SPLIT_FOURIERTOR = nloctyp.lt.0
634 nloctyp = iabs(nloctyp)
636 read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
637 read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
638 itype2loc(ntyp1)=nloctyp
639 iloctyp(nloctyp)=ntyp1
641 itype2loc(-i)=-itype2loc(i)
650 iloctyp(-i)=-iloctyp(i)
652 c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
653 c write (iout,*) "nloctyp",nloctyp,
654 c & " iloctyp",(iloctyp(i),i=0,nloctyp)
657 c write (iout,*) "NEWCORR",i
658 read (ifourier,*,end=115,err=115)
661 read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
664 c write (iout,*) "NEWCORR BNEW1"
665 c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
668 read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
671 c write (iout,*) "NEWCORR BNEW2"
672 c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
674 read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
675 read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
677 c write (iout,*) "NEWCORR CCNEW"
678 c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
680 read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
681 read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
683 c write (iout,*) "NEWCORR DDNEW"
684 c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
688 read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
692 c write (iout,*) "NEWCORR EENEW1"
693 c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
695 read (ifourier,*,end=115,err=115) e0new(ii,i)
697 c write (iout,*) (e0new(ii,i),ii=1,3)
699 c write (iout,*) "NEWCORR EENEW"
702 ccnew(ii,1,i)=ccnew(ii,1,i)/2
703 ccnew(ii,2,i)=ccnew(ii,2,i)/2
704 ddnew(ii,1,i)=ddnew(ii,1,i)/2
705 ddnew(ii,2,i)=ddnew(ii,2,i)/2
710 bnew1(ii,1,-i)= bnew1(ii,1,i)
711 bnew1(ii,2,-i)=-bnew1(ii,2,i)
712 bnew2(ii,1,-i)= bnew2(ii,1,i)
713 bnew2(ii,2,-i)=-bnew2(ii,2,i)
716 c ccnew(ii,1,i)=ccnew(ii,1,i)/2
717 c ccnew(ii,2,i)=ccnew(ii,2,i)/2
718 c ddnew(ii,1,i)=ddnew(ii,1,i)/2
719 c ddnew(ii,2,i)=ddnew(ii,2,i)/2
720 ccnew(ii,1,-i)=ccnew(ii,1,i)
721 ccnew(ii,2,-i)=-ccnew(ii,2,i)
722 ddnew(ii,1,-i)=ddnew(ii,1,i)
723 ddnew(ii,2,-i)=-ddnew(ii,2,i)
725 e0new(1,-i)= e0new(1,i)
726 e0new(2,-i)=-e0new(2,i)
727 e0new(3,-i)=-e0new(3,i)
729 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
730 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
731 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
732 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
736 write (iout,'(a)') "Coefficients of the multibody terms"
737 do i=-nloctyp+1,nloctyp-1
738 write (iout,*) "Type: ",onelet(iloctyp(i))
739 write (iout,*) "Coefficients of the expansion of B1"
741 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
743 write (iout,*) "Coefficients of the expansion of B2"
745 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
747 write (iout,*) "Coefficients of the expansion of C"
748 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
749 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
750 write (iout,*) "Coefficients of the expansion of D"
751 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
752 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
753 write (iout,*) "Coefficients of the expansion of E"
754 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
757 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
762 IF (SPLIT_FOURIERTOR) THEN
764 c write (iout,*) "NEWCORR TOR",i
765 read (ifourier,*,end=115,err=115)
768 read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i)
771 c write (iout,*) "NEWCORR BNEW1 TOR"
772 c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
775 read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i)
778 c write (iout,*) "NEWCORR BNEW2 TOR"
779 c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
781 read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i)
782 read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i)
784 c write (iout,*) "NEWCORR CCNEW TOR"
785 c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
787 read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i)
788 read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i)
790 c write (iout,*) "NEWCORR DDNEW TOR"
791 c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
795 read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i)
799 c write (iout,*) "NEWCORR EENEW1 TOR"
800 c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
802 read (ifourier,*,end=115,err=115) e0newtor(ii,i)
804 c write (iout,*) (e0newtor(ii,i),ii=1,3)
806 c write (iout,*) "NEWCORR EENEW TOR"
809 ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
810 ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
811 ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
812 ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
817 bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
818 bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
819 bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
820 bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
823 ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
824 ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
825 ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
826 ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
828 e0newtor(1,-i)= e0newtor(1,i)
829 e0newtor(2,-i)=-e0newtor(2,i)
830 e0newtor(3,-i)=-e0newtor(3,i)
832 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
833 eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
834 eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
835 eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
840 & "Single-body coefficients of the torsional potentials"
841 do i=-nloctyp+1,nloctyp-1
842 write (iout,*) "Type: ",onelet(iloctyp(i))
843 write (iout,*) "Coefficients of the expansion of B1tor"
845 write (iout,'(3hB1(,i1,1h),3f10.5)')
846 & j,(bnew1tor(k,j,i),k=1,3)
848 write (iout,*) "Coefficients of the expansion of B2tor"
850 write (iout,'(3hB2(,i1,1h),3f10.5)')
851 & j,(bnew2tor(k,j,i),k=1,3)
853 write (iout,*) "Coefficients of the expansion of Ctor"
854 write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
855 write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
856 write (iout,*) "Coefficients of the expansion of Dtor"
857 write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
858 write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
859 write (iout,*) "Coefficients of the expansion of Etor"
860 write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
863 write (iout,'(1hE,2i1,2f10.5)')
864 & j,k,(eenewtor(l,j,k,i),l=1,2)
870 do i=-nloctyp+1,nloctyp-1
873 bnew1tor(ii,j,i)=bnew1(ii,j,i)
878 bnew2tor(ii,j,i)=bnew2(ii,j,i)
882 ccnewtor(ii,1,i)=ccnew(ii,1,i)
883 ccnewtor(ii,2,i)=ccnew(ii,2,i)
884 ddnewtor(ii,1,i)=ddnew(ii,1,i)
885 ddnewtor(ii,2,i)=ddnew(ii,2,i)
888 ENDIF !SPLIT_FOURIER_TOR
891 & write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
893 read (ifourier,*,end=115,err=115)
894 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
896 write (iout,*) 'Type ',onelet(iloctyp(i))
897 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
911 c B1tilde(1,i) = b(3)
912 c B1tilde(2,i) =-b(5)
913 c B1tilde(1,-i) =-b(3)
914 c B1tilde(2,-i) =b(5)
921 cc B1tilde(1,i) = b(3,i)
922 cc B1tilde(2,i) =-b(5,i)
923 C B1tilde(1,-i) =-b(3,i)
924 C B1tilde(2,-i) =b(5,i)
925 cc b1tilde(1,i)=0.0d0
926 cc b1tilde(2,i)=0.0d0
938 CCold(1,1,-i)= b(7,i)
939 CCold(2,2,-i)=-b(7,i)
940 CCold(2,1,-i)=-b(9,i)
941 CCold(1,2,-i)=-b(9,i)
946 c Ctilde(1,1,i)= CCold(1,1,i)
947 c Ctilde(1,2,i)= CCold(1,2,i)
948 c Ctilde(2,1,i)=-CCold(2,1,i)
949 c Ctilde(2,2,i)=-CCold(2,2,i)
951 c Ctilde(1,1,i)=0.0d0
952 c Ctilde(1,2,i)=0.0d0
953 c Ctilde(2,1,i)=0.0d0
954 c Ctilde(2,2,i)=0.0d0
959 DDold(1,1,-i)= b(6,i)
960 DDold(2,2,-i)=-b(6,i)
961 DDold(2,1,-i)=-b(8,i)
962 DDold(1,2,-i)=-b(8,i)
967 c Dtilde(1,1,i)= DD(1,1,i)
968 c Dtilde(1,2,i)= DD(1,2,i)
969 c Dtilde(2,1,i)=-DD(2,1,i)
970 c Dtilde(2,2,i)=-DD(2,2,i)
972 c Dtilde(1,1,i)=0.0d0
973 c Dtilde(1,2,i)=0.0d0
974 c Dtilde(2,1,i)=0.0d0
975 c Dtilde(2,2,i)=0.0d0
976 EEold(1,1,i)= b(10,i)+b(11,i)
977 EEold(2,2,i)=-b(10,i)+b(11,i)
978 EEold(2,1,i)= b(12,i)-b(13,i)
979 EEold(1,2,i)= b(12,i)+b(13,i)
980 EEold(1,1,-i)= b(10,i)+b(11,i)
981 EEold(2,2,-i)=-b(10,i)+b(11,i)
982 EEold(2,1,-i)=-b(12,i)+b(13,i)
983 EEold(1,2,-i)=-b(12,i)-b(13,i)
984 write(iout,*) "TU DOCHODZE"
990 c ee(2,1,i)=ee(1,2,i)
995 &"Coefficients of the cumulants (independent of valence angles)"
996 do i=-nloctyp+1,nloctyp-1
997 write (iout,*) 'Type ',onelet(iloctyp(i))
999 write(iout,'(2f10.5)') B(3,i),B(5,i)
1001 write(iout,'(2f10.5)') B(2,i),B(4,i)
1004 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
1008 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
1012 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1020 C Read torsional parameters in old format
1022 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1023 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1024 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1029 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1035 write (iout,'(/a/)') 'Torsional constants:'
1038 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1039 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1045 C Read torsional parameters
1047 IF (TOR_MODE.eq.0) THEN
1049 read (itorp,*,end=113,err=113) ntortyp
1050 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1052 itype2loc(i)=itortyp(i)
1055 itype2loc(-i)=-itype2loc(i)
1057 itortyp(ntyp1)=ntortyp
1060 itortyp(i)=-itortyp(-i)
1062 write (iout,*) 'ntortyp',ntortyp
1064 do j=-ntortyp+1,ntortyp-1
1065 read (itorp,*,end=113,err=113) nterm(i,j,iblock),
1067 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1068 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1071 do k=1,nterm(i,j,iblock)
1072 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
1074 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1075 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1076 v0ij=v0ij+si*v1(k,i,j,iblock)
1078 c write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
1079 c write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
1080 c &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
1082 do k=1,nlor(i,j,iblock)
1083 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
1084 & vlor2(k,i,j),vlor3(k,i,j)
1085 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1088 v0(-i,-j,iblock)=v0ij
1094 write (iout,'(/a/)') 'Parameters of the SCCOR potentials:'
1097 do j=-ntortyp+1,ntortyp-1
1098 write (iout,*) 'ityp',i,' jtyp',j
1099 write (iout,*) 'Fourier constants'
1100 do k=1,nterm(i,j,iblock)
1101 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
1104 write (iout,*) 'Lorenz constants'
1105 do k=1,nlor(i,j,iblock)
1106 write (iout,'(3(1pe15.5))')
1107 & vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1115 C 6/23/01 Read parameters for double torsionals
1119 do j=-ntortyp+1,ntortyp-1
1120 do k=-ntortyp+1,ntortyp-1
1121 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1122 c write (iout,*) "OK onelett",
1125 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
1126 & .or. t3.ne.toronelet(k)) then
1127 write (iout,*) "Error in double torsional parameter file",
1130 call MPI_Finalize(Ierror)
1132 stop "Error in double torsional parameter file"
1134 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
1135 & ntermd_2(i,j,k,iblock)
1136 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1137 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1138 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
1139 & ntermd_1(i,j,k,iblock))
1140 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
1141 & ntermd_1(i,j,k,iblock))
1142 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
1143 & ntermd_1(i,j,k,iblock))
1144 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
1145 & ntermd_1(i,j,k,iblock))
1146 C Martix of D parameters for one dimesional foureir series
1147 do l=1,ntermd_1(i,j,k,iblock)
1148 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1149 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1150 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1151 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1152 c write(iout,*) "whcodze" ,
1153 c & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1155 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
1156 & v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
1157 & v2s(m,l,i,j,k,iblock),
1158 & m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1159 C Martix of D parameters for two dimesional fourier series
1160 do l=1,ntermd_2(i,j,k,iblock)
1162 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1163 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1164 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1165 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1174 write (iout,*) 'Constants for double torsionals'
1177 do j=-ntortyp+1,ntortyp-1
1178 do k=-ntortyp+1,ntortyp-1
1179 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
1180 & ' nsingle',ntermd_1(i,j,k,iblock),
1181 & ' ndouble',ntermd_2(i,j,k,iblock)
1183 write (iout,*) 'Single angles:'
1184 do l=1,ntermd_1(i,j,k,iblock)
1185 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
1186 & v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
1187 & v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
1188 & v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1191 write (iout,*) 'Pairs of angles:'
1192 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1193 do l=1,ntermd_2(i,j,k,iblock)
1194 write (iout,'(i5,20f10.5)')
1195 & l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1198 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1199 do l=1,ntermd_2(i,j,k,iblock)
1200 write (iout,'(i5,20f10.5)')
1201 & l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
1202 & (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1211 ELSE IF (TOR_MODE.eq.1) THEN
1213 C read valence-torsional parameters
1214 read (itorp,*,end=121,err=121) ntortyp
1216 write (iout,*) "Valence-torsional parameters read in ntortyp",
1218 read (itorp,*,end=121,err=121) (itortyp(i),i=1,ntyp)
1219 write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
1222 itype2loc(i)=itortyp(i)
1226 itortyp(i)=-itortyp(-i)
1228 do i=-ntortyp+1,ntortyp-1
1229 do j=-ntortyp+1,ntortyp-1
1230 C first we read the cos and sin gamma parameters
1231 read (itorp,'(13x,a)',end=121,err=121) string
1232 write (iout,*) i,j,string
1233 read (itorp,*,end=121,err=121)
1234 & nterm_kcc(j,i),nterm_kcc_Tb(j,i)
1235 C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
1236 do k=1,nterm_kcc(j,i)
1237 do l=1,nterm_kcc_Tb(j,i)
1238 do ll=1,nterm_kcc_Tb(j,i)
1239 read (itorp,*,end=121,err=121) ii,jj,kk,
1240 & v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1248 c AL 4/8/16: Calculate coefficients from one-body parameters
1251 itortyp(i)=itype2loc(i)
1254 &"Val-tor parameters calculated from cumulant coefficients ntortyp"
1256 do i=-ntortyp+1,ntortyp-1
1257 do j=-ntortyp+1,ntortyp-1
1260 do k=1,nterm_kcc_Tb(j,i)
1261 do l=1,nterm_kcc_Tb(j,i)
1262 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)
1263 & +bnew1tor(k,2,i)*bnew2tor(l,2,j)
1264 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)
1265 & +bnew1tor(k,2,i)*bnew2tor(l,1,j)
1268 do k=1,nterm_kcc_Tb(j,i)
1269 do l=1,nterm_kcc_Tb(j,i)
1271 v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1272 & -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1273 v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1274 & +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1276 v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1277 & -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1278 v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1279 & +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1283 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)
1287 write (iout,*) "TOR_MODE>1 only with NEWCORR"
1292 if (tor_mode.gt.0 .and. lprint) then
1293 c Print valence-torsional parameters
1295 & "Parameters of the valence-torsional potentials"
1296 do i=-ntortyp+1,ntortyp-1
1297 do j=-ntortyp+1,ntortyp-1
1298 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
1299 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
1300 do k=1,nterm_kcc(j,i)
1301 do l=1,nterm_kcc_Tb(j,i)
1302 do ll=1,nterm_kcc_Tb(j,i)
1303 write (iout,'(3i5,2f15.4)')
1304 & k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1313 C Read of Side-chain backbone correlation parameters
1314 C Modified 11 May 2012 by Adasko
1317 read (isccor,*,end=119,err=119) nsccortyp
1319 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1321 isccortyp(i)=-isccortyp(-i)
1323 iscprol=isccortyp(20)
1324 c write (iout,*) 'ntortyp',ntortyp
1326 cc maxinter is maximum interaction sites
1330 read (isccor,*,end=119,err=119)
1331 &nterm_sccor(i,j),nlor_sccor(i,j)
1337 nterm_sccor(-i,j)=nterm_sccor(i,j)
1338 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1339 nterm_sccor(i,-j)=nterm_sccor(i,j)
1340 do k=1,nterm_sccor(i,j)
1341 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
1343 if (j.eq.iscprol) then
1344 if (i.eq.isccortyp(10)) then
1345 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1346 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1348 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
1349 & +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1350 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
1351 & +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1352 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1353 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1354 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1355 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1358 if (i.eq.isccortyp(10)) then
1359 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1360 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1362 if (j.eq.isccortyp(10)) then
1363 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1364 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1366 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1367 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1368 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1369 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1370 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1371 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1375 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1376 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1377 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1378 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1381 do k=1,nlor_sccor(i,j)
1382 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
1383 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1384 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
1385 &(1+vlor3sccor(k,i,j)**2)
1387 v0sccor(l,i,j)=v0ijsccor
1388 v0sccor(l,-i,j)=v0ijsccor1
1389 v0sccor(l,i,-j)=v0ijsccor2
1390 v0sccor(l,-i,-j)=v0ijsccor3
1396 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1397 c write (iout,*) 'ntortyp',ntortyp
1399 cc maxinter is maximum interaction sites
1403 read (isccor,*,end=119,err=119)
1404 & nterm_sccor(i,j),nlor_sccor(i,j)
1408 do k=1,nterm_sccor(i,j)
1409 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
1411 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1414 do k=1,nlor_sccor(i,j)
1415 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
1416 & vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1417 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
1418 &(1+vlor3sccor(k,i,j)**2)
1420 v0sccor(l,i,j)=v0ijsccor
1428 write (iout,'(/a/)') 'Torsional constants:'
1432 write (iout,*) 'ityp',i,' jtyp',j
1433 write (iout,*) 'Fourier constants'
1434 do k=1,nterm_sccor(i,j)
1435 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1437 write (iout,*) 'Lorenz constants'
1438 do k=1,nlor_sccor(i,j)
1439 write (iout,'(3(1pe15.5))')
1440 & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1448 C Read electrostatic-interaction parameters
1452 write (iout,'(/a)') 'Electrostatic interaction constants:'
1453 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1454 & 'IT','JT','APP','BPP','AEL6','AEL3'
1456 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1457 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1458 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1459 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1464 app (i,j)=epp(i,j)*rri*rri
1465 bpp (i,j)=-2.0D0*epp(i,j)*rri
1466 ael6(i,j)=elpp6(i,j)*4.2D0**6
1467 ael3(i,j)=elpp3(i,j)*4.2D0**3
1469 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1470 & ael6(i,j),ael3(i,j)
1475 C Read side-chain interaction parameters.
1477 read (isidep,*,end=117,err=117) ipot,expon
1478 if (ipot.lt.1 .or. ipot.gt.5) then
1479 write (iout,'(2a)') 'Error while reading SC interaction',
1480 & 'potential file - unknown potential type.'
1482 call MPI_Finalize(Ierror)
1488 & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1489 & ', exponents are ',expon,2*expon
1490 goto (10,20,30,30,40) ipot
1491 C----------------------- LJ potential ---------------------------------
1492 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1493 & (sigma0(i),i=1,ntyp)
1495 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1496 write (iout,'(a/)') 'The epsilon array:'
1497 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1498 write (iout,'(/a)') 'One-body parameters:'
1499 write (iout,'(a,4x,a)') 'residue','sigma'
1500 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1503 C----------------------- LJK potential --------------------------------
1504 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1505 & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1507 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1508 write (iout,'(a/)') 'The epsilon array:'
1509 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1510 write (iout,'(/a)') 'One-body parameters:'
1511 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1512 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1516 C---------------------- GB or BP potential -----------------------------
1518 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1520 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
1521 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
1522 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
1523 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
1524 C now we start reading lipid
1526 read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
1527 if (lprint) write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp)
1529 C print *,"WARNING!!"
1531 C epslip(i,j)=epslip(i,j)+0.05d0
1534 if (lprint) write(iout,*) epslip(1,1),"OK?"
1535 C For the GB potential convert sigma'**2 into chi'
1538 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1542 write (iout,'(/a/)') 'Parameters of the BP potential:'
1543 write (iout,'(a/)') 'The epsilon array:'
1544 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1545 write (iout,'(/a)') 'One-body parameters:'
1546 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',
1548 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1549 & chip(i),alp(i),i=1,ntyp)
1552 C--------------------- GBV potential -----------------------------------
1553 c 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1554 c & (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1555 c & (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1557 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1559 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
1560 read (isidep,*,end=117,err=117)(rr0(i),i=1,ntyp)
1561 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
1562 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
1563 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
1565 read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
1566 if (lprint) write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp)
1569 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1572 write (iout,'(/a/)') 'Parameters of the GBV potential:'
1573 write (iout,'(a/)') 'The epsilon array:'
1574 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1575 write (iout,'(/a)') 'One-body parameters:'
1576 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',
1577 & 's||/s_|_^2',' chip ',' alph '
1578 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1579 & sigii(i),chip(i),alp(i),i=1,ntyp)
1581 C now we start reading lipid
1584 C-----------------------------------------------------------------------
1585 C Calculate the "working" parameters of SC interactions.
1589 epslip(i,j)=epslip(j,i)
1594 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1595 sigma(j,i)=sigma(i,j)
1596 rs0(i,j)=dwa16*sigma(i,j)
1600 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1601 & 'Working parameters of the SC interactions:',
1602 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1607 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1616 sigeps=dsign(1.0D0,epsij)
1618 aa_aq(i,j)=epsij*rrij*rrij
1619 bb_aq(i,j)=-sigeps*epsij*rrij
1620 aa_aq(j,i)=aa_aq(i,j)
1621 bb_aq(j,i)=bb_aq(i,j)
1622 epsijlip=epslip(i,j)
1623 sigeps=dsign(1.0D0,epsijlip)
1624 epsijlip=dabs(epsijlip)
1625 aa_lip(i,j)=epsijlip*rrij*rrij
1626 bb_lip(i,j)=-sigeps*epsijlip*rrij
1627 aa_lip(j,i)=aa_lip(i,j)
1628 bb_lip(j,i)=bb_lip(i,j)
1630 sigt1sq=sigma0(i)**2
1631 sigt2sq=sigma0(j)**2
1634 ratsig1=sigt2sq/sigt1sq
1635 ratsig2=1.0D0/ratsig1
1636 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1637 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1638 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1642 c if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1643 sigmaii(i,j)=rsum_max
1644 sigmaii(j,i)=rsum_max
1646 c sigmaii(i,j)=r0(i,j)
1647 c sigmaii(j,i)=r0(i,j)
1649 cd write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1650 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1651 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1652 augm(i,j)=epsij*r_augm**(2*expon)
1653 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1660 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1661 & restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1662 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1667 write(iout,*) "tube param"
1668 read(itube,*) epspeptube,sigmapeptube
1669 sigmapeptube=sigmapeptube**6
1670 sigeps=dsign(1.0D0,epspeptube)
1671 epspeptube=dabs(epspeptube)
1672 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
1673 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
1674 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
1676 read(itube,*) epssctube,sigmasctube
1677 sigmasctube=sigmasctube**6
1678 sigeps=dsign(1.0D0,epssctube)
1679 epssctube=dabs(epssctube)
1680 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
1681 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
1682 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
1688 C Define the SC-p interaction constants (hard-coded; old style)
1691 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1693 c aad(i,1)=0.3D0*4.0D0**12
1694 C Following line for constants currently implemented
1695 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1696 aad(i,1)=1.5D0*4.0D0**12
1697 c aad(i,1)=0.17D0*5.6D0**12
1699 C "Soft" SC-p repulsion
1701 C Following line for constants currently implemented
1702 c aad(i,1)=0.3D0*4.0D0**6
1703 C "Hard" SC-p repulsion
1704 bad(i,1)=3.0D0*4.0D0**6
1705 c bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1714 C 8/9/01 Read the SC-p interaction constants from file
1717 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1720 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1721 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1722 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1723 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1727 write (iout,*) "Parameters of SC-p interactions:"
1729 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1730 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1736 C Define the constants of the disulfide bridge
1740 c Old arbitrary potential - commented out.
1745 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1746 c energy surface of diethyl disulfide.
1747 c A. Liwo and U. Kozlowska, 11/24/03
1763 C if(me.eq.king) then
1764 C write (iout,'(/a)') "Disulfide bridge parameters:"
1765 C write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1766 C write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1767 C write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1768 C write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1771 if (shield_mode.gt.0) then
1772 C VSolvSphere the volume of solving sphere
1773 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1774 C there will be no distinction between proline peptide group and normal peptide
1775 C group in case of shielding parameters
1776 c write (iout,*) "rpp(1,1)",rpp(1,1)," pi",pi
1777 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1778 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1779 write (iout,*) "VSolvSphere",VSolvSphere,"VSolvSphere_div",
1781 C long axis of side chain
1783 long_r_sidechain(i)=vbldsc0(1,i)
1784 short_r_sidechain(i)=sigma0(i)
1789 111 write (iout,*) "Error reading bending energy parameters."
1791 112 write (iout,*) "Error reading rotamer energy parameters."
1793 113 write (iout,*) "Error reading torsional energy parameters."
1795 114 write (iout,*) "Error reading double torsional energy parameters."
1798 & "Error reading cumulant (multibody energy) parameters."
1800 116 write (iout,*) "Error reading electrostatic energy parameters."
1802 1161 write (iout,*) "Error reading electrostatic energy parameters.Lip"
1804 117 write (iout,*) "Error reading side chain interaction parameters."
1806 118 write (iout,*) "Error reading SCp interaction parameters."
1808 119 write (iout,*) "Error reading SCCOR parameters"
1810 121 write (iout,*) "Error in Czybyshev parameters"
1813 call MPI_Finalize(Ierror)
1820 subroutine getenv_loc(var, val)
1821 character(*) var, val
1824 character(2000) line
1827 open (196,file='env',status='old',readonly,shared)
1829 c write(*,*)'looking for ',var
1830 10 read(196,*,err=11,end=11)line
1831 iread=index(line,var)
1832 c write(*,*)iread,' ',var,' ',line
1833 if (iread.eq.0) go to 10
1834 c write(*,*)'---> ',line
1840 iread=iread+ilen(var)+1
1841 read (line(iread:),*,err=12,end=12) val
1842 c write(*,*)'OK: ',var,' = ',val
1848 #elif (defined CRAY)
1849 integer lennam,lenval,ierror
1851 c getenv using a POSIX call, useful on the T3D
1852 c Sept 1996, comment out error check on advice of H. Pritchard
1855 if(lennam.le.0) stop '--error calling getenv--'
1856 call pxfgetenv(var,lennam,val,lenval,ierror)
1857 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1859 call getenv(var,val)
1861 C set the variables used for shielding effect
1862 C if (shield_mode.gt.0) then
1863 C VSolvSphere the volume of solving sphere
1865 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1866 C there will be no distinction between proline peptide group and normal peptide
1867 C group in case of shielding parameters
1868 C VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1869 C VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1870 C long axis of side chain
1872 C long_r_sidechain(i)=vbldsc0(1,i)
1873 C short_r_sidechain(i)=sigma0(i)
1875 C lets set the buffor value