wham and cluster_wham Adam's new constr_dist multichain
[unres.git] / source / unres / src_MD-M / parmread.F
1       subroutine parmread
2 C
3 C Read the parameters of the probability distributions of the virtual-bond
4 C valence angles and the side chains and energy parameters.
5 C
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
8 C parameters.
9 C
10       implicit real*8 (a-h,o-z)
11       include 'DIMENSIONS'
12 #ifdef MPI
13       include "mpif.h"
14       integer IERROR
15 #endif
16       include 'COMMON.IOUNITS'
17       include 'COMMON.CHAIN'
18       include 'COMMON.INTERACT'
19       include 'COMMON.GEO'
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'
27       include 'COMMON.MD'
28       include 'COMMON.SETUP'
29       character*1 t1,t2,t3
30       character*1 onelett(4) /"G","A","P","D"/
31       character*1 toronelet(-2:2) /"p","a","G","A","P"/
32       logical lprint,LaTeX
33       dimension blower(3,3,maxlob)
34 C      dimension b(13)
35       character*3 lancuch,ucase
36 C
37 C For printing parameters after they are read set the following in the UNRES
38 C C-shell script:
39 C
40 C setenv PRINT_PARM YES
41 C
42 C To print parameters in LaTeX format rather than as ASCII tables:
43 C
44 C setenv LATEX YES
45 C
46       call getenv_loc("PRINT_PARM",lancuch)
47       lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
48       call getenv_loc("LATEX",lancuch)
49       LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
50 C
51       dwa16=2.0d0**(1.0d0/6.0d0)
52       itypro=20
53 C Assign virtual-bond length
54       vbl=3.8D0
55       vblinv=1.0D0/vbl
56       vblinv2=vblinv*vblinv
57 c
58 c Read the virtual-bond parameters, masses, and moments of inertia
59 c and Stokes' radii of the peptide group and side chains
60 c
61 #ifdef CRYST_BOND
62       read (ibond,*) vbldp0,vbldpdum,akp,mp,ip,pstok
63       do i=1,ntyp
64         nbondterm(i)=1
65         read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
66         dsc(i) = vbldsc0(1,i)
67         if (i.eq.10) then
68           dsc_inv(i)=0.0D0
69         else
70           dsc_inv(i)=1.0D0/dsc(i)
71         endif
72       enddo
73 #else
74       read (ibond,*) junk,vbldp0,vbldpdum,akp,rjunk,mp,ip,pstok
75       do i=1,ntyp
76         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
77      &   j=1,nbondterm(i)),msc(i),isc(i),restok(i)
78         dsc(i) = vbldsc0(1,i)
79         if (i.eq.10) then
80           dsc_inv(i)=0.0D0
81         else
82           dsc_inv(i)=1.0D0/dsc(i)
83         endif
84       enddo
85 #endif
86       if (lprint) then
87         write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
88         write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',
89      &   'inertia','Pstok'
90         write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
91         do i=1,ntyp
92           write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
93      &      vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
94           do j=2,nbondterm(i)
95             write (iout,'(13x,3f10.5)')
96      &        vbldsc0(j,i),aksc(j,i),abond0(j,i)
97           enddo
98         enddo
99       endif
100 C reading lipid parameters
101        read(iliptranpar,*) pepliptran
102        do i=1,ntyp
103        read(iliptranpar,*) liptranene(i)
104        enddo
105        close(iliptranpar)
106 #ifdef CRYST_THETA
107 C
108 C Read the parameters of the probability distribution/energy expression 
109 C of the virtual-bond valence angles theta
110 C
111       do i=1,ntyp
112         read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),
113      &    (bthet(j,i,1,1),j=1,2)
114         read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
115         read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
116         read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
117         sigc0(i)=sigc0(i)**2
118       enddo
119       do i=1,ntyp
120       athet(1,i,1,-1)=athet(1,i,1,1)
121       athet(2,i,1,-1)=athet(2,i,1,1)
122       bthet(1,i,1,-1)=-bthet(1,i,1,1)
123       bthet(2,i,1,-1)=-bthet(2,i,1,1)
124       athet(1,i,-1,1)=-athet(1,i,1,1)
125       athet(2,i,-1,1)=-athet(2,i,1,1)
126       bthet(1,i,-1,1)=bthet(1,i,1,1)
127       bthet(2,i,-1,1)=bthet(2,i,1,1)
128       enddo
129       do i=-ntyp,-1
130       a0thet(i)=a0thet(-i)
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)
135       athet(1,i,-1,1)=athet(1,-i,1,1)
136       athet(2,i,-1,1)=-athet(2,-i,1,1)
137       bthet(1,i,-1,1)=-bthet(1,-i,1,1)
138       bthet(2,i,-1,1)=bthet(2,-i,1,1)
139       athet(1,i,1,-1)=-athet(1,-i,1,1)
140       athet(2,i,1,-1)=athet(2,-i,1,1)
141       bthet(1,i,1,-1)=bthet(1,-i,1,1)
142       bthet(2,i,1,-1)=-bthet(2,-i,1,1)
143       theta0(i)=theta0(-i)
144       sig0(i)=sig0(-i)
145       sigc0(i)=sigc0(-i)
146        do j=0,3
147         polthet(j,i)=polthet(j,-i)
148        enddo
149        do j=1,3
150          gthet(j,i)=gthet(j,-i)
151        enddo
152       enddo
153
154       close (ithep)
155       if (lprint) then
156       if (.not.LaTeX) then
157         write (iout,'(a)') 
158      &    'Parameters of the virtual-bond valence angles:'
159         write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
160      & '    ATHETA0   ','         A1   ','        A2    ',
161      & '        B1    ','         B2   '        
162         do i=1,ntyp
163           write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
164      &        a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
165         enddo
166         write (iout,'(/a/9x,5a/79(1h-))') 
167      & 'Parameters of the expression for sigma(theta_c):',
168      & '     ALPH0    ','      ALPH1   ','     ALPH2    ',
169      & '     ALPH3    ','    SIGMA0C   '        
170         do i=1,ntyp
171           write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
172      &      (polthet(j,i),j=0,3),sigc0(i) 
173         enddo
174         write (iout,'(/a/9x,5a/79(1h-))') 
175      & 'Parameters of the second gaussian:',
176      & '    THETA0    ','     SIGMA0   ','        G1    ',
177      & '        G2    ','         G3   '        
178         do i=1,ntyp
179           write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
180      &       sig0(i),(gthet(j,i),j=1,3)
181         enddo
182        else
183         write (iout,'(a)') 
184      &    'Parameters of the virtual-bond valence angles:'
185         write (iout,'(/a/9x,5a/79(1h-))') 
186      & 'Coefficients of expansion',
187      & '     theta0   ','    a1*10^2   ','   a2*10^2    ',
188      & '   b1*10^1    ','    b2*10^1   '        
189         do i=1,ntyp
190           write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
191      &        a0thet(i),(100*athet(j,i,1,1),j=1,2),
192      &        (10*bthet(j,i,1,1),j=1,2)
193         enddo
194         write (iout,'(/a/9x,5a/79(1h-))') 
195      & 'Parameters of the expression for sigma(theta_c):',
196      & ' alpha0       ','  alph1       ',' alph2        ',
197      & ' alhp3        ','   sigma0c    '        
198         do i=1,ntyp
199           write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
200      &      (polthet(j,i),j=0,3),sigc0(i) 
201         enddo
202         write (iout,'(/a/9x,5a/79(1h-))') 
203      & 'Parameters of the second gaussian:',
204      & '    theta0    ','  sigma0*10^2 ','      G1*10^-1',
205      & '        G2    ','   G3*10^1    '        
206         do i=1,ntyp
207           write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
208      &       100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
209         enddo
210       endif
211       endif
212 #else 
213
214 C Read the parameters of Utheta determined from ab initio surfaces
215 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
216 C
217       read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,
218      &  ntheterm3,nsingle,ndouble
219       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
220       read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
221       do i=-ntyp1,-1
222         ithetyp(i)=-ithetyp(-i)
223       enddo
224       do iblock=1,2
225       do i=-maxthetyp,maxthetyp
226         do j=-maxthetyp,maxthetyp
227           do k=-maxthetyp,maxthetyp
228             aa0thet(i,j,k,iblock)=0.0d0
229             do l=1,ntheterm
230               aathet(l,i,j,k,iblock)=0.0d0
231             enddo
232             do l=1,ntheterm2
233               do m=1,nsingle
234                 bbthet(m,l,i,j,k,iblock)=0.0d0
235                 ccthet(m,l,i,j,k,iblock)=0.0d0
236                 ddthet(m,l,i,j,k,iblock)=0.0d0
237                 eethet(m,l,i,j,k,iblock)=0.0d0
238               enddo
239             enddo
240             do l=1,ntheterm3
241               do m=1,ndouble
242                 do mm=1,ndouble
243                  ffthet(mm,m,l,i,j,k,iblock)=0.0d0
244                  ggthet(mm,m,l,i,j,k,iblock)=0.0d0
245                 enddo
246               enddo
247             enddo
248           enddo
249         enddo
250       enddo
251       enddo
252 c VAR:iblock means terminally blocking group 1=non-proline 2=proline
253       do iblock=1,2 
254 c VAR:ntethtyp is type of theta potentials type currently 0=glycine 
255 c VAR:1=non-glicyne non-proline 2=proline
256 c VAR:negative values for D-aminoacid
257       do i=0,nthetyp
258         do j=-nthetyp,nthetyp
259           do k=-nthetyp,nthetyp
260             read (ithep,'(6a)',end=111,err=111) res1
261             read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
262 c VAR: aa0thet is variable describing the average value of Foureir
263 c VAR: expansion series
264 c VAR: aathet is foureir expansion in theta/2 angle for full formula
265 c VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
266 Condens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
267             read (ithep,*,end=111,err=111)
268      &(aathet(l,i,j,k,iblock),l=1,ntheterm)
269             read (ithep,*,end=111,err=111)
270      &       ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
271      &        (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
272      &        (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
273      &        (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),
274      &        ll=1,ntheterm2)
275             read (ithep,*,end=111,err=111)
276      &      (((ffthet(llll,lll,ll,i,j,k,iblock),
277      &         ffthet(lll,llll,ll,i,j,k,iblock),
278      &         ggthet(llll,lll,ll,i,j,k,iblock),
279      &         ggthet(lll,llll,ll,i,j,k,iblock),
280      &         llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
281           enddo
282         enddo
283       enddo
284 C
285 C For dummy ends assign glycine-type coefficients of theta-only terms; the
286 C coefficients of theta-and-gamma-dependent terms are zero.
287 C IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
288 C RECOMENTDED AFTER VERSION 3.3)
289 c      do i=1,nthetyp
290 c        do j=1,nthetyp
291 c          do l=1,ntheterm
292 c            aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
293 c            aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
294 c          enddo
295 c          aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
296 c          aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
297 c        enddo
298 c        do l=1,ntheterm
299 c          aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
300 c        enddo
301 c        aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
302 c      enddo
303 c      enddo
304 C AND COMMENT THE LOOPS BELOW
305       do i=1,nthetyp
306         do j=1,nthetyp
307           do l=1,ntheterm
308             aathet(l,i,j,nthetyp+1,iblock)=0.0d0
309             aathet(l,nthetyp+1,i,j,iblock)=0.0d0
310           enddo
311           aa0thet(i,j,nthetyp+1,iblock)=0.0d0
312           aa0thet(nthetyp+1,i,j,iblock)=0.0d0
313         enddo
314         do l=1,ntheterm
315           aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
316         enddo
317         aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
318       enddo
319       enddo
320 C TILL HERE
321 C Substitution for D aminoacids from symmetry.
322       do iblock=1,2
323       do i=-nthetyp,0
324         do j=-nthetyp,nthetyp
325           do k=-nthetyp,nthetyp
326            aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
327            do l=1,ntheterm
328            aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock) 
329            enddo
330            do ll=1,ntheterm2
331             do lll=1,nsingle
332             bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
333             ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
334             ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
335             eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
336             enddo
337           enddo
338           do ll=1,ntheterm3
339            do lll=2,ndouble
340             do llll=1,lll-1
341             ffthet(llll,lll,ll,i,j,k,iblock)=
342      &      ffthet(llll,lll,ll,-i,-j,-k,iblock) 
343             ffthet(lll,llll,ll,i,j,k,iblock)=
344      &      ffthet(lll,llll,ll,-i,-j,-k,iblock)
345             ggthet(llll,lll,ll,i,j,k,iblock)=
346      &      -ggthet(llll,lll,ll,-i,-j,-k,iblock)
347             ggthet(lll,llll,ll,i,j,k,iblock)=
348      &      -ggthet(lll,llll,ll,-i,-j,-k,iblock)      
349             enddo !ll
350            enddo  !lll  
351           enddo   !llll
352          enddo    !k
353         enddo     !j
354        enddo      !i
355       enddo       !iblock
356 C
357 C Control printout of the coefficients of virtual-bond-angle potentials
358 C
359       if (lprint) then
360         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
361         do iblock=1,2
362         do i=1,nthetyp+1
363           do j=1,nthetyp+1
364             do k=1,nthetyp+1
365               write (iout,'(//4a)') 
366      &         'Type ',onelett(i),onelett(j),onelett(k) 
367               write (iout,'(//a,10x,a)') " l","a[l]"
368               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
369               write (iout,'(i2,1pe15.5)')
370      &           (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
371             do l=1,ntheterm2
372               write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') 
373      &          "b",l,"c",l,"d",l,"e",l
374               do m=1,nsingle
375                 write (iout,'(i2,4(1pe15.5))') m,
376      &          bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
377      &          ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
378               enddo
379             enddo
380             do l=1,ntheterm3
381               write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
382      &          "f+",l,"f-",l,"g+",l,"g-",l
383               do m=2,ndouble
384                 do n=1,m-1
385                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
386      &              ffthet(n,m,l,i,j,k,iblock),
387      &              ffthet(m,n,l,i,j,k,iblock),
388      &              ggthet(n,m,l,i,j,k,iblock),
389      &              ggthet(m,n,l,i,j,k,iblock)
390                 enddo
391               enddo
392             enddo
393           enddo
394         enddo
395       enddo
396       enddo
397       call flush(iout)
398       endif
399 c      write (2,*) "Start reading THETA_PDB",ithep_pdb
400       do i=1,ntyp
401 c      write (2,*) 'i=',i
402         read (ithep_pdb,*,err=111,end=111)
403      &     a0thet(i),(athet(j,i,1,1),j=1,2),
404      &    (bthet(j,i,1,1),j=1,2)
405         read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
406         read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
407         read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
408         sigc0(i)=sigc0(i)**2
409       enddo
410       do i=1,ntyp
411       athet(1,i,1,-1)=athet(1,i,1,1)
412       athet(2,i,1,-1)=athet(2,i,1,1)
413       bthet(1,i,1,-1)=-bthet(1,i,1,1)
414       bthet(2,i,1,-1)=-bthet(2,i,1,1)
415       athet(1,i,-1,1)=-athet(1,i,1,1)
416       athet(2,i,-1,1)=-athet(2,i,1,1)
417       bthet(1,i,-1,1)=bthet(1,i,1,1)
418       bthet(2,i,-1,1)=bthet(2,i,1,1)
419       enddo
420       do i=-ntyp,-1
421       a0thet(i)=a0thet(-i)
422       athet(1,i,-1,-1)=athet(1,-i,1,1)
423       athet(2,i,-1,-1)=-athet(2,-i,1,1)
424       bthet(1,i,-1,-1)=bthet(1,-i,1,1)
425       bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
426       athet(1,i,-1,1)=athet(1,-i,1,1)
427       athet(2,i,-1,1)=-athet(2,-i,1,1)
428       bthet(1,i,-1,1)=-bthet(1,-i,1,1)
429       bthet(2,i,-1,1)=bthet(2,-i,1,1)
430       athet(1,i,1,-1)=-athet(1,-i,1,1)
431       athet(2,i,1,-1)=athet(2,-i,1,1)
432       bthet(1,i,1,-1)=bthet(1,-i,1,1)
433       bthet(2,i,1,-1)=-bthet(2,-i,1,1)
434       theta0(i)=theta0(-i)
435       sig0(i)=sig0(-i)
436       sigc0(i)=sigc0(-i)
437        do j=0,3
438         polthet(j,i)=polthet(j,-i)
439        enddo
440        do j=1,3
441          gthet(j,i)=gthet(j,-i)
442        enddo
443       enddo
444 c      write (2,*) "End reading THETA_PDB"
445       close (ithep_pdb)
446 #endif
447       close(ithep)
448 #ifdef CRYST_SC
449 C
450 C Read the parameters of the probability distribution/energy expression
451 C of the side chains.
452 C
453       do i=1,ntyp
454         read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
455         if (i.eq.10) then
456           dsc_inv(i)=0.0D0
457         else
458           dsc_inv(i)=1.0D0/dsc(i)
459         endif
460         if (i.ne.10) then
461         do j=1,nlob(i)
462           do k=1,3
463             do l=1,3
464               blower(l,k,j)=0.0D0
465             enddo
466           enddo
467         enddo  
468         bsc(1,i)=0.0D0
469         read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
470      &    ((blower(k,l,1),l=1,k),k=1,3)
471         censc(1,1,-i)=censc(1,1,i)
472         censc(2,1,-i)=censc(2,1,i)
473         censc(3,1,-i)=-censc(3,1,i)
474         do j=2,nlob(i)
475           read (irotam,*,end=112,err=112) bsc(j,i)
476           read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
477      &                                 ((blower(k,l,j),l=1,k),k=1,3)
478         censc(1,j,-i)=censc(1,j,i)
479         censc(2,j,-i)=censc(2,j,i)
480         censc(3,j,-i)=-censc(3,j,i)
481 C BSC is amplitude of Gaussian
482         enddo
483         do j=1,nlob(i)
484           do k=1,3
485             do l=1,k
486               akl=0.0D0
487               do m=1,3
488                 akl=akl+blower(k,m,j)*blower(l,m,j)
489               enddo
490               gaussc(k,l,j,i)=akl
491               gaussc(l,k,j,i)=akl
492              if (((k.eq.3).and.(l.ne.3))
493      &        .or.((l.eq.3).and.(k.ne.3))) then
494                 gaussc(k,l,j,-i)=-akl
495                 gaussc(l,k,j,-i)=-akl
496               else
497                 gaussc(k,l,j,-i)=akl
498                 gaussc(l,k,j,-i)=akl
499               endif
500             enddo
501           enddo 
502         enddo
503         endif
504       enddo
505       close (irotam)
506       if (lprint) then
507         write (iout,'(/a)') 'Parameters of side-chain local geometry'
508         do i=1,ntyp
509           nlobi=nlob(i)
510           if (nlobi.gt.0) then
511             if (LaTeX) then
512               write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
513      &         ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
514                write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
515      &                             'log h',(bsc(j,i),j=1,nlobi)
516                write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
517      &        'x',((censc(k,j,i),k=1,3),j=1,nlobi)
518               do k=1,3
519                 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
520      &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
521               enddo
522             else
523               write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
524               write (iout,'(a,f10.4,4(16x,f10.4))')
525      &                             'Center  ',(bsc(j,i),j=1,nlobi)
526               write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
527      &           j=1,nlobi)
528               write (iout,'(a)')
529             endif
530           endif
531         enddo
532       endif
533 #else
534
535 C Read scrot parameters for potentials determined from all-atom AM1 calculations
536 C added by Urszula Kozlowska 07/11/2007
537 C
538       do i=1,ntyp
539         read (irotam,*,end=112,err=112) 
540        if (i.eq.10) then 
541          read (irotam,*,end=112,err=112) 
542        else
543          do j=1,65
544            read(irotam,*,end=112,err=112) sc_parmin(j,i)
545          enddo  
546        endif
547       enddo
548 C
549 C Read the parameters of the probability distribution/energy expression
550 C of the side chains.
551 C
552 c      write (2,*) "Start reading ROTAM_PDB"
553       do i=1,ntyp
554         read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
555         if (i.eq.10) then
556           dsc_inv(i)=0.0D0
557         else
558           dsc_inv(i)=1.0D0/dsc(i)
559         endif
560         if (i.ne.10) then
561         do j=1,nlob(i)
562           do k=1,3
563             do l=1,3
564               blower(l,k,j)=0.0D0
565             enddo
566           enddo
567         enddo
568         bsc(1,i)=0.0D0
569         read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
570      &    ((blower(k,l,1),l=1,k),k=1,3)
571         do j=2,nlob(i)
572           read (irotam_pdb,*,end=112,err=112) bsc(j,i)
573           read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
574      &                                 ((blower(k,l,j),l=1,k),k=1,3)
575         enddo
576         do j=1,nlob(i)
577           do k=1,3
578             do l=1,k
579               akl=0.0D0
580               do m=1,3
581                 akl=akl+blower(k,m,j)*blower(l,m,j)
582               enddo
583               gaussc(k,l,j,i)=akl
584               gaussc(l,k,j,i)=akl
585             enddo
586           enddo
587         enddo
588         endif
589       enddo
590       close (irotam_pdb)
591 c      write (2,*) "End reading ROTAM_PDB"
592 #endif
593       close(irotam)
594
595 #ifdef CRYST_TOR
596 C
597 C Read torsional parameters in old format
598 C
599       read (itorp,*,end=113,err=113) ntortyp,nterm_old
600       if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
601       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
602       do i=1,ntortyp
603         do j=1,ntortyp
604           read (itorp,'(a)')
605           do k=1,nterm_old
606             read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i) 
607           enddo
608         enddo
609       enddo
610       close (itorp)
611       if (lprint) then
612         write (iout,'(/a/)') 'Torsional constants:'
613         do i=1,ntortyp
614           do j=1,ntortyp
615             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
616             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
617           enddo
618         enddo
619       endif
620 #else
621 C
622 C Read torsional parameters
623 C
624       read (itorp,*,end=113,err=113) ntortyp
625       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
626       do iblock=1,2
627       do i=-ntyp,-1
628        itortyp(i)=-itortyp(-i)
629       enddo
630 c      write (iout,*) 'ntortyp',ntortyp
631       do i=0,ntortyp-1
632         do j=-ntortyp+1,ntortyp-1
633           read (itorp,*,end=113,err=113) nterm(i,j,iblock),
634      &          nlor(i,j,iblock)
635           nterm(-i,-j,iblock)=nterm(i,j,iblock)
636           nlor(-i,-j,iblock)=nlor(i,j,iblock)
637           v0ij=0.0d0
638           si=-1.0d0
639           do k=1,nterm(i,j,iblock)
640             read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
641      &      v2(k,i,j,iblock)
642             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
643             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
644             v0ij=v0ij+si*v1(k,i,j,iblock)
645             si=-si
646 c         write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
647 c         write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
648 c      &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
649           enddo
650           do k=1,nlor(i,j,iblock)
651             read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
652      &        vlor2(k,i,j),vlor3(k,i,j)
653             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
654           enddo
655           v0(i,j,iblock)=v0ij
656           v0(-i,-j,iblock)=v0ij
657         enddo
658       enddo
659       enddo
660       close (itorp)
661       if (lprint) then
662         write (iout,'(/a/)') 'Torsional constants:'
663         do iblock=1,2
664         do i=1,ntortyp
665           do j=1,ntortyp
666             write (iout,*) 'ityp',i,' jtyp',j
667             write (iout,*) 'Fourier constants'
668             do k=1,nterm(i,j,iblock)
669               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
670      &        v2(k,i,j,iblock)
671             enddo
672             write (iout,*) 'Lorenz constants'
673             do k=1,nlor(i,j,iblock)
674               write (iout,'(3(1pe15.5))')
675      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
676             enddo
677           enddo
678         enddo
679         enddo
680       endif
681
682 C
683 C 6/23/01 Read parameters for double torsionals
684 C
685       do iblock=1,2
686       do i=0,ntortyp-1
687         do j=-ntortyp+1,ntortyp-1
688           do k=-ntortyp+1,ntortyp-1
689             read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
690 c              write (iout,*) "OK onelett",
691 c     &         i,j,k,t1,t2,t3
692
693             if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
694      &        .or. t3.ne.toronelet(k)) then
695               write (iout,*) "Error in double torsional parameter file",
696      &         i,j,k,t1,t2,t3
697 #ifdef MPI
698               call MPI_Finalize(Ierror)
699 #endif
700                stop "Error in double torsional parameter file"
701             endif
702            read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
703      &         ntermd_2(i,j,k,iblock)
704             ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
705             ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
706             read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
707      &         ntermd_1(i,j,k,iblock))
708             read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
709      &         ntermd_1(i,j,k,iblock))
710             read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
711      &         ntermd_1(i,j,k,iblock))
712             read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
713      &         ntermd_1(i,j,k,iblock))
714 C Martix of D parameters for one dimesional foureir series
715             do l=1,ntermd_1(i,j,k,iblock)
716              v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
717              v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
718              v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
719              v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
720 c            write(iout,*) "whcodze" ,
721 c     & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
722             enddo
723             read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
724      &         v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
725      &         v2s(m,l,i,j,k,iblock),
726      &         m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
727 C Martix of D parameters for two dimesional fourier series
728             do l=1,ntermd_2(i,j,k,iblock)
729              do m=1,l-1
730              v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
731              v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
732              v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
733              v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
734              enddo!m
735             enddo!l
736           enddo!k
737         enddo!j
738       enddo!i
739       enddo!iblock
740       if (lprint) then
741       write (iout,*)
742       write (iout,*) 'Constants for double torsionals'
743       do iblock=1,2
744       do i=0,ntortyp-1
745         do j=-ntortyp+1,ntortyp-1
746           do k=-ntortyp+1,ntortyp-1
747             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
748      &        ' nsingle',ntermd_1(i,j,k,iblock),
749      &        ' ndouble',ntermd_2(i,j,k,iblock)
750             write (iout,*)
751             write (iout,*) 'Single angles:'
752             do l=1,ntermd_1(i,j,k,iblock)
753               write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
754      &           v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
755      &           v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
756      &           v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
757             enddo
758             write (iout,*)
759             write (iout,*) 'Pairs of angles:'
760             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
761             do l=1,ntermd_2(i,j,k,iblock)
762               write (iout,'(i5,20f10.5)')
763      &         l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
764             enddo
765             write (iout,*)
766             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
767             do l=1,ntermd_2(i,j,k,iblock)
768               write (iout,'(i5,20f10.5)')
769      &         l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
770      &         (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
771             enddo
772             write (iout,*)
773           enddo
774         enddo
775       enddo
776       enddo
777       endif
778 #endif
779 C Read of Side-chain backbone correlation parameters
780 C Modified 11 May 2012 by Adasko
781 CCC
782 C
783       read (isccor,*,end=119,err=119) nsccortyp
784 #ifdef SCCORPDB
785       read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
786       do i=-ntyp,-1
787         isccortyp(i)=-isccortyp(-i)
788       enddo
789       iscprol=isccortyp(20)
790 c      write (iout,*) 'ntortyp',ntortyp
791       maxinter=3
792 cc maxinter is maximum interaction sites
793       do l=1,maxinter
794       do i=1,nsccortyp
795         do j=1,nsccortyp
796           read (isccor,*,end=119,err=119)
797      &nterm_sccor(i,j),nlor_sccor(i,j)
798           v0ijsccor=0.0d0
799           v0ijsccor1=0.0d0
800           v0ijsccor2=0.0d0
801           v0ijsccor3=0.0d0
802           si=-1.0d0
803           nterm_sccor(-i,j)=nterm_sccor(i,j)
804           nterm_sccor(-i,-j)=nterm_sccor(i,j)
805           nterm_sccor(i,-j)=nterm_sccor(i,j)
806           do k=1,nterm_sccor(i,j)
807             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
808      &    ,v2sccor(k,l,i,j)
809             if (j.eq.iscprol) then
810              if (i.eq.isccortyp(10)) then
811              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
812              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
813              else
814              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
815      &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
816              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
817      &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
818              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
819              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
820              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
821              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
822              endif
823             else
824              if (i.eq.isccortyp(10)) then
825              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
826              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
827              else
828                if (j.eq.isccortyp(10)) then
829              v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
830              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
831                else
832              v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
833              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
834              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
835              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
836              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
837              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
838                 endif
839                endif
840             endif
841             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
842             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
843             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
844             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
845             si=-si
846           enddo
847           do k=1,nlor_sccor(i,j)
848             read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
849      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
850             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
851      &(1+vlor3sccor(k,i,j)**2)
852           enddo
853           v0sccor(l,i,j)=v0ijsccor
854           v0sccor(l,-i,j)=v0ijsccor1
855           v0sccor(l,i,-j)=v0ijsccor2
856           v0sccor(l,-i,-j)=v0ijsccor3         
857         enddo
858       enddo
859       enddo
860       close (isccor)
861 #else
862       read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
863 c      write (iout,*) 'ntortyp',ntortyp
864       maxinter=3
865 cc maxinter is maximum interaction sites
866       do l=1,maxinter
867       do i=1,nsccortyp
868         do j=1,nsccortyp
869           read (isccor,*,end=119,err=119)
870      & nterm_sccor(i,j),nlor_sccor(i,j)
871           v0ijsccor=0.0d0
872           si=-1.0d0
873
874           do k=1,nterm_sccor(i,j)
875             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
876      &    ,v2sccor(k,l,i,j)
877             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
878             si=-si
879           enddo
880           do k=1,nlor_sccor(i,j)
881             read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
882      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
883             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
884      &(1+vlor3sccor(k,i,j)**2)
885           enddo
886            v0sccor(l,i,j)=v0ijsccor
887         enddo
888       enddo
889       enddo
890       close (isccor)
891
892 #endif      
893       if (lprint) then
894         write (iout,'(/a/)') 'Torsional constants:'
895         do l=1,maxinter
896         do i=1,nsccortyp
897           do j=1,nsccortyp
898             write (iout,*) 'ityp',i,' jtyp',j
899             write (iout,*) 'Fourier constants'
900             do k=1,nterm_sccor(i,j)
901       write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
902             enddo
903             write (iout,*) 'Lorenz constants'
904             do k=1,nlor_sccor(i,j)
905               write (iout,'(3(1pe15.5))')
906      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
907             enddo
908           enddo
909         enddo
910         enddo
911       endif
912
913 C
914 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
915 C         interaction energy of the Gly, Ala, and Pro prototypes.
916 C
917       if (lprint) then
918         write (iout,*)
919         write (iout,*) "Coefficients of the cumulants"
920       endif
921       read (ifourier,*) nloctyp
922
923       do i=0,nloctyp-1
924         read (ifourier,*,end=115,err=115)
925         read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
926 #ifdef NEWCORR
927         read (ifourier,*,end=115,err=115) (bnew1(ii,1,i),ii=1,3)
928         read (ifourier,*,end=115,err=115) (bnew2(ii,1,i),ii=1,3)
929         read (ifourier,*,end=115,err=115) (bnew1(ii,2,i),ii=1,1)
930         read (ifourier,*,end=115,err=115) (bnew2(ii,2,i),ii=1,1)
931         read (ifourier,*,end=115,err=115) (eenew(ii,i),ii=1,1)
932 #endif 
933         if (lprint) then
934         write (iout,*) 'Type',i
935         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
936         endif
937         b(3,-i)= b(3,i)
938         b(5,-i)=-b(5,i)
939         b(4,-i)=-b(4,i)
940         b(2,-i)= b(2,i)     
941         B1(1,i)  = b(3,i)
942         B1(2,i)  = b(5,i)
943         B1(1,-i) = b(3,i)
944         B1(2,-i) = -b(5,i)
945 c        b1(1,i)=0.0d0
946 c        b1(2,i)=0.0d0
947         B1tilde(1,i) = b(3,i)
948         B1tilde(2,i) =-b(5,i)
949         B1tilde(1,-i) =-b(3,i)
950         B1tilde(2,-i) =b(5,i)
951 c        b1tilde(1,i)=0.0d0
952 c        b1tilde(2,i)=0.0d0
953         B2(1,i)  = b(2,i)
954         B2(2,i)  = b(4,i)
955         B2(1,-i)  =b(2,i)
956         B2(2,-i)  =-b(4,i)
957 c        b2(1,i)=0.0d0
958 c        b2(2,i)=0.0d0
959         CC(1,1,i)= b(7,i)
960         CC(2,2,i)=-b(7,i)
961         CC(2,1,i)= b(9,i)
962         CC(1,2,i)= b(9,i)
963         CC(1,1,-i)= b(7,i)
964         CC(2,2,-i)=-b(7,i)
965         CC(2,1,-i)=-b(9,i)
966         CC(1,2,-i)=-b(9,i)
967 c        CC(1,1,i)=0.0d0
968 c        CC(2,2,i)=0.0d0
969 c        CC(2,1,i)=0.0d0
970 c        CC(1,2,i)=0.0d0
971         Ctilde(1,1,i)=b(7,i)
972         Ctilde(1,2,i)=b(9,i)
973         Ctilde(2,1,i)=-b(9,i)
974         Ctilde(2,2,i)=b(7,i)
975         Ctilde(1,1,-i)=b(7,i)
976         Ctilde(1,2,-i)=-b(9,i)
977         Ctilde(2,1,-i)=b(9,i)
978         Ctilde(2,2,-i)=b(7,i)
979
980 c        Ctilde(1,1,i)=0.0d0
981 c        Ctilde(1,2,i)=0.0d0
982 c        Ctilde(2,1,i)=0.0d0
983 c        Ctilde(2,2,i)=0.0d0
984         DD(1,1,i)= b(6,i)
985         DD(2,2,i)=-b(6,i)
986         DD(2,1,i)= b(8,i)
987         DD(1,2,i)= b(8,i)
988         DD(1,1,-i)= b(6,i)
989         DD(2,2,-i)=-b(6,i)
990         DD(2,1,-i)=-b(8,i)
991         DD(1,2,-i)=-b(8,i)
992 c        DD(1,1,i)=0.0d0
993 c        DD(2,2,i)=0.0d0
994 c        DD(2,1,i)=0.0d0
995 c        DD(1,2,i)=0.0d0
996         Dtilde(1,1,i)=b(6,i)
997         Dtilde(1,2,i)=b(8,i)
998         Dtilde(2,1,i)=-b(8,i)
999         Dtilde(2,2,i)=b(6,i)
1000         Dtilde(1,1,-i)=b(6,i)
1001         Dtilde(1,2,-i)=-b(8,i)
1002         Dtilde(2,1,-i)=b(8,i)
1003         Dtilde(2,2,-i)=b(6,i)
1004
1005 c        Dtilde(1,1,i)=0.0d0
1006 c        Dtilde(1,2,i)=0.0d0
1007 c        Dtilde(2,1,i)=0.0d0
1008 c        Dtilde(2,2,i)=0.0d0
1009         EEold(1,1,i)= b(10,i)+b(11,i)
1010         EEold(2,2,i)=-b(10,i)+b(11,i)
1011         EEold(2,1,i)= b(12,i)-b(13,i)
1012         EEold(1,2,i)= b(12,i)+b(13,i)
1013         EEold(1,1,-i)= b(10,i)+b(11,i)
1014         EEold(2,2,-i)=-b(10,i)+b(11,i)
1015         EEold(2,1,-i)=-b(12,i)+b(13,i)
1016         EEold(1,2,-i)=-b(12,i)-b(13,i)
1017 c        ee(1,1,i)=1.0d0
1018 c        ee(2,2,i)=1.0d0
1019 c        ee(2,1,i)=0.0d0
1020 c        ee(1,2,i)=0.0d0
1021 c        ee(2,1,i)=ee(1,2,i)
1022       enddo
1023 c      lprint=.true.
1024       if (lprint) then
1025       do i=1,nloctyp
1026         write (iout,*) 'Type',i
1027         write (iout,*) 'B1'
1028         write(iout,*) B1(1,i),B1(2,i)
1029         write (iout,*) 'B2'
1030         write(iout,*) B2(1,i),B2(2,i)
1031         write (iout,*) 'CC'
1032         do j=1,2
1033           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1034         enddo
1035         write(iout,*) 'DD'
1036         do j=1,2
1037           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1038         enddo
1039         write(iout,*) 'EE'
1040         do j=1,2
1041           write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1042         enddo
1043       enddo
1044       endif
1045 c      lprint=.false.
1046
1047
1048 C Read electrostatic-interaction parameters
1049 C
1050       if (lprint) then
1051         write (iout,*)
1052         write (iout,'(/a)') 'Electrostatic interaction constants:'
1053         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') 
1054      &            'IT','JT','APP','BPP','AEL6','AEL3'
1055       endif
1056       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1057       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1058       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1059       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1060       close (ielep)
1061       do i=1,2
1062         do j=1,2
1063         rri=rpp(i,j)**6
1064         app (i,j)=epp(i,j)*rri*rri 
1065         bpp (i,j)=-2.0D0*epp(i,j)*rri
1066         ael6(i,j)=elpp6(i,j)*4.2D0**6
1067         ael3(i,j)=elpp3(i,j)*4.2D0**3
1068 c        lprint=.true.
1069         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1070      &                    ael6(i,j),ael3(i,j)
1071 c        lprint=.false.
1072         enddo
1073       enddo
1074 C
1075 C Read side-chain interaction parameters.
1076 C
1077       read (isidep,*,end=117,err=117) ipot,expon
1078       if (ipot.lt.1 .or. ipot.gt.5) then
1079         write (iout,'(2a)') 'Error while reading SC interaction',
1080      &               'potential file - unknown potential type.'
1081 #ifdef MPI
1082         call MPI_Finalize(Ierror)
1083 #endif
1084         stop
1085       endif
1086       expon2=expon/2
1087       if(me.eq.king)
1088      & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1089      & ', exponents are ',expon,2*expon 
1090       goto (10,20,30,30,40) ipot
1091 C----------------------- LJ potential ---------------------------------
1092    10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1093      &   (sigma0(i),i=1,ntyp)
1094       if (lprint) then
1095         write (iout,'(/a/)') 'Parameters of the LJ potential:'
1096         write (iout,'(a/)') 'The epsilon array:'
1097         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1098         write (iout,'(/a)') 'One-body parameters:'
1099         write (iout,'(a,4x,a)') 'residue','sigma'
1100         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1101       endif
1102       goto 50
1103 C----------------------- LJK potential --------------------------------
1104    20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1105      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1106       if (lprint) then
1107         write (iout,'(/a/)') 'Parameters of the LJK potential:'
1108         write (iout,'(a/)') 'The epsilon array:'
1109         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1110         write (iout,'(/a)') 'One-body parameters:'
1111         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
1112         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1113      &        i=1,ntyp)
1114       endif
1115       goto 50
1116 C---------------------- GB or BP potential -----------------------------
1117    30 do i=1,ntyp
1118        read (isidep,*,end=116,err=116)(eps(i,j),j=i,ntyp)
1119       enddo
1120       read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
1121       read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
1122       read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
1123       read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
1124 C now we start reading lipid
1125       do i=1,ntyp
1126        read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
1127 C       print *,"WARNING!!"
1128 C       do j=1,ntyp
1129 C       epslip(i,j)=epslip(i,j)+0.05d0
1130 C       enddo
1131       enddo
1132 c      write(iout,*) epslip(1,1),"OK?"
1133 C For the GB potential convert sigma'**2 into chi'
1134       if (ipot.eq.4) then
1135         do i=1,ntyp
1136           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1137         enddo
1138       endif
1139       if (lprint) then
1140         write (iout,'(/a/)') 'Parameters of the BP potential:'
1141         write (iout,'(a/)') 'The epsilon array:'
1142         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1143         write (iout,'(/a)') 'One-body parameters:'
1144         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
1145      &       '    chip  ','    alph  '
1146         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1147      &                     chip(i),alp(i),i=1,ntyp)
1148       endif
1149       goto 50
1150 C--------------------- GBV potential -----------------------------------
1151    40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1152      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1153      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1154       if (lprint) then
1155         write (iout,'(/a/)') 'Parameters of the GBV potential:'
1156         write (iout,'(a/)') 'The epsilon array:'
1157         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1158         write (iout,'(/a)') 'One-body parameters:'
1159         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
1160      &      's||/s_|_^2','    chip  ','    alph  '
1161         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1162      &           sigii(i),chip(i),alp(i),i=1,ntyp)
1163       endif
1164    50 continue
1165       close (isidep)
1166 C-----------------------------------------------------------------------
1167 C Calculate the "working" parameters of SC interactions.
1168       do i=2,ntyp
1169         do j=1,i-1
1170           eps(i,j)=eps(j,i)
1171           epslip(i,j)=epslip(j,i)
1172         enddo
1173       enddo
1174       do i=1,ntyp
1175         do j=i,ntyp
1176           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1177           sigma(j,i)=sigma(i,j)
1178           rs0(i,j)=dwa16*sigma(i,j)
1179           rs0(j,i)=rs0(i,j)
1180         enddo
1181       enddo
1182       if (lprint) write (iout,'(/a/10x,7a/72(1h-))') 
1183      & 'Working parameters of the SC interactions:',
1184      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
1185      & '  chi1   ','   chi2   ' 
1186       do i=1,ntyp
1187         do j=i,ntyp
1188           epsij=eps(i,j)
1189           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1190             rrij=sigma(i,j)
1191           else
1192             rrij=rr0(i)+rr0(j)
1193           endif
1194           r0(i,j)=rrij
1195           r0(j,i)=rrij
1196           rrij=rrij**expon
1197           epsij=eps(i,j)
1198           sigeps=dsign(1.0D0,epsij)
1199           epsij=dabs(epsij)
1200           aa_aq(i,j)=epsij*rrij*rrij
1201           bb_aq(i,j)=-sigeps*epsij*rrij
1202           aa_aq(j,i)=aa_aq(i,j)
1203           bb_aq(j,i)=bb_aq(i,j)
1204           epsijlip=epslip(i,j)
1205           sigeps=dsign(1.0D0,epsijlip)
1206           epsijlip=dabs(epsijlip)
1207           aa_lip(i,j)=epsijlip*rrij*rrij
1208           bb_lip(i,j)=-sigeps*epsijlip*rrij
1209           aa_lip(j,i)=aa_lip(i,j)
1210           bb_lip(j,i)=bb_lip(i,j)
1211           if (ipot.gt.2) then
1212             sigt1sq=sigma0(i)**2
1213             sigt2sq=sigma0(j)**2
1214             sigii1=sigii(i)
1215             sigii2=sigii(j)
1216             ratsig1=sigt2sq/sigt1sq
1217             ratsig2=1.0D0/ratsig1
1218             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1219             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1220             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1221           else
1222             rsum_max=sigma(i,j)
1223           endif
1224 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1225             sigmaii(i,j)=rsum_max
1226             sigmaii(j,i)=rsum_max 
1227 c         else
1228 c           sigmaii(i,j)=r0(i,j)
1229 c           sigmaii(j,i)=r0(i,j)
1230 c         endif
1231 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1232           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1233             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1234             augm(i,j)=epsij*r_augm**(2*expon)
1235 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1236             augm(j,i)=augm(i,j)
1237           else
1238             augm(i,j)=0.0D0
1239             augm(j,i)=0.0D0
1240           endif
1241           if (lprint) then
1242             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') 
1243      &      restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1244      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1245           endif
1246         enddo
1247       enddo
1248 #ifdef OLDSCP
1249 C
1250 C Define the SC-p interaction constants (hard-coded; old style)
1251 C
1252       do i=1,ntyp
1253 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1254 C helix formation)
1255 c       aad(i,1)=0.3D0*4.0D0**12
1256 C Following line for constants currently implemented
1257 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1258         aad(i,1)=1.5D0*4.0D0**12
1259 c       aad(i,1)=0.17D0*5.6D0**12
1260         aad(i,2)=aad(i,1)
1261 C "Soft" SC-p repulsion
1262         bad(i,1)=0.0D0
1263 C Following line for constants currently implemented
1264 c       aad(i,1)=0.3D0*4.0D0**6
1265 C "Hard" SC-p repulsion
1266         bad(i,1)=3.0D0*4.0D0**6
1267 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1268         bad(i,2)=bad(i,1)
1269 c       aad(i,1)=0.0D0
1270 c       aad(i,2)=0.0D0
1271 c       bad(i,1)=1228.8D0
1272 c       bad(i,2)=1228.8D0
1273       enddo
1274 #else
1275 C
1276 C 8/9/01 Read the SC-p interaction constants from file
1277 C
1278       do i=1,ntyp
1279         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1280       enddo
1281       do i=1,ntyp
1282         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1283         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1284         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1285         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1286       enddo
1287 c      lprint=.true.
1288       if (lprint) then
1289         write (iout,*) "Parameters of SC-p interactions:"
1290         do i=1,ntyp
1291           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1292      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1293         enddo
1294       endif
1295 c      lprint=.false.
1296 #endif
1297 C
1298 C Define the constants of the disulfide bridge
1299 C
1300       ebr=-5.50D0
1301 c
1302 c Old arbitrary potential - commented out.
1303 c
1304 c      dbr= 4.20D0
1305 c      fbr= 3.30D0
1306 c
1307 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1308 c energy surface of diethyl disulfide.
1309 c A. Liwo and U. Kozlowska, 11/24/03
1310 c
1311       D0CM = 3.78d0
1312       AKCM = 15.1d0
1313       AKTH = 11.0d0
1314       AKCT = 12.0d0
1315       V1SS =-1.08d0
1316       V2SS = 7.61d0
1317       V3SS = 13.7d0
1318 c      akcm=0.0d0
1319 c      akth=0.0d0
1320 c      akct=0.0d0
1321 c      v1ss=0.0d0
1322 c      v2ss=0.0d0
1323 c      v3ss=0.0d0
1324       
1325       if(me.eq.king) then
1326       write (iout,'(/a)') "Disulfide bridge parameters:"
1327       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1328       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1329       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1330       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1331      &  ' v3ss:',v3ss
1332       endif
1333       return
1334   111 write (iout,*) "Error reading bending energy parameters."
1335       goto 999
1336   112 write (iout,*) "Error reading rotamer energy parameters."
1337       goto 999
1338   113 write (iout,*) "Error reading torsional energy parameters."
1339       goto 999
1340   114 write (iout,*) "Error reading double torsional energy parameters."
1341       goto 999
1342   115 write (iout,*) 
1343      &  "Error reading cumulant (multibody energy) parameters."
1344       goto 999
1345   116 write (iout,*) "Error reading electrostatic energy parameters."
1346       goto 999
1347  1161 write (iout,*) "Error reading electrostatic energy parameters.Lip"
1348       goto 999
1349   117 write (iout,*) "Error reading side chain interaction parameters."
1350       goto 999
1351   118 write (iout,*) "Error reading SCp interaction parameters."
1352       goto 999
1353   119 write (iout,*) "Error reading SCCOR parameters"
1354   999 continue
1355 #ifdef MPI
1356       call MPI_Finalize(Ierror)
1357 #endif
1358       stop
1359       return
1360       end
1361
1362
1363       subroutine getenv_loc(var, val)
1364       character(*) var, val
1365
1366 #ifdef WINIFL
1367       character(2000) line
1368       external ilen
1369
1370       open (196,file='env',status='old',readonly,shared)
1371       iread=0
1372 c      write(*,*)'looking for ',var
1373 10    read(196,*,err=11,end=11)line
1374       iread=index(line,var)
1375 c      write(*,*)iread,' ',var,' ',line
1376       if (iread.eq.0) go to 10 
1377 c      write(*,*)'---> ',line
1378 11    continue
1379       if(iread.eq.0) then
1380 c       write(*,*)'CHUJ'
1381        val=''
1382       else
1383        iread=iread+ilen(var)+1
1384        read (line(iread:),*,err=12,end=12) val
1385 c       write(*,*)'OK: ',var,' = ',val
1386       endif
1387       close(196)
1388       return
1389 12    val=''
1390       close(196)
1391 #elif (defined CRAY)
1392       integer lennam,lenval,ierror
1393 c
1394 c        getenv using a POSIX call, useful on the T3D
1395 c        Sept 1996, comment out error check on advice of H. Pritchard
1396 c
1397       lennam = len(var)
1398       if(lennam.le.0) stop '--error calling getenv--'
1399       call pxfgetenv(var,lennam,val,lenval,ierror)
1400 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1401 #else
1402       call getenv(var,val)
1403 #endif
1404
1405       return
1406       end