Adam's changes
[unres.git] / source / wham / src-M-Pawel / parmread.F
1       subroutine parmread(iparm,*)
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       implicit real*8 (a-h,o-z)
7       include 'DIMENSIONS'
8       include 'DIMENSIONS.ZSCOPT'
9       include 'DIMENSIONS.FREE'
10       include 'COMMON.IOUNITS'
11       include 'COMMON.CHAIN'
12       include 'COMMON.INTERACT'
13       include 'COMMON.GEO'
14       include 'COMMON.LOCAL'
15       include 'COMMON.TORSION'
16       include 'COMMON.FFIELD'
17       include 'COMMON.NAMES'
18       include 'COMMON.SBRIDGE'
19       include 'COMMON.WEIGHTS'
20       include 'COMMON.ENEPS'
21       include 'COMMON.SCCOR'
22       include 'COMMON.SCROT'
23       include 'COMMON.FREE'
24       character*1 t1,t2,t3
25       character*1 onelett(4) /"G","A","P","D"/
26       character*1 toronelet(-2:2) /"p","a","G","A","P"/
27       logical lprint
28       dimension blower(3,3,maxlob)
29       character*800 controlcard
30       character*256 bondname_t,thetname_t,rotname_t,torname_t,
31      & tordname_t,fouriername_t,elename_t,sidename_t,scpname_t,
32      & sccorname_t
33       integer ilen
34       external ilen
35       character*16 key
36       integer iparm
37       double precision ip,mp
38 C
39 C Body
40 C
41 C Set LPRINT=.TRUE. for debugging
42       dwa16=2.0d0**(1.0d0/6.0d0)
43       lprint=.false.
44       itypro=20
45 C Assign virtual-bond length
46       vbl=3.8D0
47       vblinv=1.0D0/vbl
48       vblinv2=vblinv*vblinv
49       call card_concat(controlcard,.true.)
50       wname(4)="WCORRH"
51       do i=1,n_ene
52         key = wname(i)(:ilen(wname(i)))
53         call reada(controlcard,key(:ilen(key)),ww(i),1.0d0)
54       enddo
55
56       write (iout,*) "iparm",iparm," myparm",myparm
57 c If reading not own parameters, skip assignment
58
59       if (iparm.eq.myparm .or. .not.separate_parset) then
60
61 c
62 c Setup weights for UNRES
63 c
64       wsc=ww(1)
65       wscp=ww(2)
66       welec=ww(3)
67       wcorr=ww(4)
68       wcorr5=ww(5)
69       wcorr6=ww(6)
70       wel_loc=ww(7)
71       wturn3=ww(8)
72       wturn4=ww(9)
73       wturn6=ww(10)
74       wang=ww(11)
75       wscloc=ww(12)
76       wtor=ww(13)
77       wtor_d=ww(14)
78       wvdwpp=ww(16)
79       wbond=ww(18)
80       wsccor=ww(19)
81
82       endif
83
84       call card_concat(controlcard,.false.)
85
86 c Return if not own parameters
87
88       if (iparm.ne.myparm .and. separate_parset) return
89
90       call reads(controlcard,"BONDPAR",bondname_t,bondname)
91       open (ibond,file=bondname_t,status='old')
92       rewind(ibond)
93       call reads(controlcard,"THETPAR",thetname_t,thetname)
94       open (ithep,file=thetname_t,status='old')
95       rewind(ithep) 
96       call reads(controlcard,"ROTPAR",rotname_t,rotname)
97       open (irotam,file=rotname_t,status='old')
98       rewind(irotam)
99       call reads(controlcard,"TORPAR",torname_t,torname)
100       open (itorp,file=torname_t,status='old')
101       rewind(itorp)
102       call reads(controlcard,"TORDPAR",tordname_t,tordname)
103       open (itordp,file=tordname_t,status='old')
104       rewind(itordp)
105       call reads(controlcard,"SCCORPAR",sccorname_t,sccorname)
106       open (isccor,file=sccorname_t,status='old')
107       rewind(isccor)
108       call reads(controlcard,"FOURIER",fouriername_t,fouriername)
109       open (ifourier,file=fouriername_t,status='old')
110       rewind(ifourier)
111       call reads(controlcard,"ELEPAR",elename_t,elename)
112       open (ielep,file=elename_t,status='old')
113       rewind(ielep)
114       call reads(controlcard,"SIDEPAR",sidename_t,sidename)
115       open (isidep,file=sidename_t,status='old')
116       rewind(isidep)
117       call reads(controlcard,"SCPPAR",scpname_t,scpname)
118       open (iscpp,file=scpname_t,status='old')
119       rewind(iscpp)
120       write (iout,*) "Parameter set:",iparm
121       write (iout,*) "Energy-term weights:"
122       do i=1,n_ene
123         write (iout,'(a16,f10.5)') wname(i),ww(i)
124       enddo
125       write (iout,*) "Sidechain potential file        : ",
126      &  sidename_t(:ilen(sidename_t))
127 #ifndef OLDSCP
128       write (iout,*) "SCp potential file              : ",
129      &  scpname_t(:ilen(scpname_t))
130 #endif  
131       write (iout,*) "Electrostatic potential file    : ",
132      &  elename_t(:ilen(elename_t))
133       write (iout,*) "Cumulant coefficient file       : ",
134      &  fouriername_t(:ilen(fouriername_t))
135       write (iout,*) "Torsional parameter file        : ",
136      &  torname_t(:ilen(torname_t))
137       write (iout,*) "Double torsional parameter file : ",
138      &  tordname_t(:ilen(tordname_t))
139       write (iout,*) "Backbone-rotamer parameter file : ",
140      &  sccorname(:ilen(sccorname))
141       write (iout,*) "Bond & inertia constant file    : ",
142      &  bondname_t(:ilen(bondname_t))
143       write (iout,*) "Bending parameter file          : ",
144      &  thetname_t(:ilen(thetname_t))
145       write (iout,*) "Rotamer parameter file          : ",
146      &  rotname_t(:ilen(rotname_t))
147
148 c
149 c Read the virtual-bond parameters, masses, and moments of inertia
150 c and Stokes' radii of the peptide group and side chains
151 c
152 #ifdef CRYST_BOND
153       read (ibond,*) vbldp0,akp
154       do i=1,ntyp
155         nbondterm(i)=1
156         read (ibond,*) vbldsc0(1,i),aksc(1,i)
157         dsc(i) = vbldsc0(1,i)
158         if (i.eq.10) then
159           dsc_inv(i)=0.0D0
160         else
161           dsc_inv(i)=1.0D0/dsc(i)
162         endif
163       enddo
164 #else
165       read (ibond,*) ijunk,vbldp0,akp,rjunk
166       do i=1,ntyp
167         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
168      &   j=1,nbondterm(i))
169         dsc(i) = vbldsc0(1,i)
170         if (i.eq.10) then
171           dsc_inv(i)=0.0D0
172         else
173           dsc_inv(i)=1.0D0/dsc(i)
174         endif
175       enddo
176 #endif
177       if (lprint) then
178         write(iout,'(/a/)')"Force constants virtual bonds:"
179         write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
180      &   'inertia','Pstok'
181         write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
182         do i=1,ntyp
183           write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
184      &      vbldsc0(1,i),aksc(1,i),abond0(1,i)
185           do j=2,nbondterm(i)
186             write (iout,'(13x,3f10.5)')
187      &        vbldsc0(j,i),aksc(j,i),abond0(j,i)
188           enddo
189         enddo
190       endif
191 #ifdef CRYST_THETA
192 C
193 C Read the parameters of the probability distribution/energy expression 
194 C of the virtual-bond valence angles theta
195 C
196       do i=1,ntyp
197         read (ithep,*) a0thet(i),(athet(j,i,1,1),j=1,2),
198      &    (bthet(j,i,1,1),j=1,2)
199         read (ithep,*) (polthet(j,i),j=0,3)
200         read (ithep,*) (gthet(j,i),j=1,3)
201         read (ithep,*) theta0(i),sig0(i),sigc0(i)
202         sigc0(i)=sigc0(i)**2
203       enddo
204       do i=1,ntyp
205       athet(1,i,1,-1)=athet(1,i,1,1)
206       athet(2,i,1,-1)=athet(2,i,1,1)
207       bthet(1,i,1,-1)=-bthet(1,i,1,1)
208       bthet(2,i,1,-1)=-bthet(2,i,1,1)
209       athet(1,i,-1,1)=-athet(1,i,1,1)
210       athet(2,i,-1,1)=-athet(2,i,1,1)
211       bthet(1,i,-1,1)=bthet(1,i,1,1)
212       bthet(2,i,-1,1)=bthet(2,i,1,1)
213       enddo
214       do i=-ntyp,-1
215       a0thet(i)=a0thet(-i)
216       athet(1,i,-1,-1)=athet(1,-i,1,1)
217       athet(2,i,-1,-1)=-athet(2,-i,1,1)
218       bthet(1,i,-1,-1)=bthet(1,-i,1,1)
219       bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
220       athet(1,i,-1,1)=athet(1,-i,1,1)
221       athet(2,i,-1,1)=-athet(2,-i,1,1)
222       bthet(1,i,-1,1)=-bthet(1,-i,1,1)
223       bthet(2,i,-1,1)=bthet(2,-i,1,1)
224       athet(1,i,1,-1)=-athet(1,-i,1,1)
225       athet(2,i,1,-1)=athet(2,-i,1,1)
226       bthet(1,i,1,-1)=bthet(1,-i,1,1)
227       bthet(2,i,1,-1)=-bthet(2,-i,1,1)
228       theta0(i)=theta0(-i)
229       sig0(i)=sig0(-i)
230       sigc0(i)=sigc0(-i)
231        do j=0,3
232         polthet(j,i)=polthet(j,-i)
233        enddo
234        do j=1,3
235          gthet(j,i)=gthet(j,-i)
236        enddo
237       enddo
238       close (ithep)
239       if (lprint) then
240 c       write (iout,'(a)') 
241 c    &    'Parameters of the virtual-bond valence angles:'
242 c       write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
243 c    & '    ATHETA0   ','         A1   ','        A2    ',
244 c    & '        B1    ','         B2   '        
245 c       do i=1,ntyp
246 c         write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
247 c    &        a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
248 c       enddo
249 c       write (iout,'(/a/9x,5a/79(1h-))') 
250 c    & 'Parameters of the expression for sigma(theta_c):',
251 c    & '     ALPH0    ','      ALPH1   ','     ALPH2    ',
252 c    & '     ALPH3    ','    SIGMA0C   '        
253 c       do i=1,ntyp
254 c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
255 c    &      (polthet(j,i),j=0,3),sigc0(i) 
256 c       enddo
257 c       write (iout,'(/a/9x,5a/79(1h-))') 
258 c    & 'Parameters of the second gaussian:',
259 c    & '    THETA0    ','     SIGMA0   ','        G1    ',
260 c    & '        G2    ','         G3   '        
261 c       do i=1,ntyp
262 c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
263 c    &       sig0(i),(gthet(j,i),j=1,3)
264 c       enddo
265         write (iout,'(a)') 
266      &    'Parameters of the virtual-bond valence angles:'
267         write (iout,'(/a/9x,5a/79(1h-))') 
268      & 'Coefficients of expansion',
269      & '     theta0   ','    a1*10^2   ','   a2*10^2    ',
270      & '   b1*10^1    ','    b2*10^1   '        
271         do i=1,ntyp
272           write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
273      &        a0thet(i),(100*athet(j,i,1,1),j=1,2),
274      &        (10*bthet(j,i,1,1),j=1,2)
275         enddo
276         write (iout,'(/a/9x,5a/79(1h-))') 
277      & 'Parameters of the expression for sigma(theta_c):',
278      & ' alpha0       ','  alph1       ',' alph2        ',
279      & ' alhp3        ','   sigma0c    '        
280         do i=1,ntyp
281           write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
282      &      (polthet(j,i),j=0,3),sigc0(i) 
283         enddo
284         write (iout,'(/a/9x,5a/79(1h-))') 
285      & 'Parameters of the second gaussian:',
286      & '    theta0    ','  sigma0*10^2 ','      G1*10^-1',
287      & '        G2    ','   G3*10^1    '        
288         do i=1,ntyp
289           write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
290      &       100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
291         enddo
292       endif
293 #else
294 C
295 C Read the parameters of Utheta determined from ab initio surfaces
296 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
297 C
298 c      write (iout,*) "tu dochodze"
299       read (ithep,*) nthetyp,ntheterm,ntheterm2,
300      &  ntheterm3,nsingle,ndouble
301       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
302       read (ithep,*) (ithetyp(i),i=1,ntyp1)
303       do i=-ntyp1,-1
304         ithetyp(i)=-ithetyp(-i)
305       enddo
306 c      write (iout,*) "tu dochodze"
307       do iblock=1,2
308       do i=-maxthetyp,maxthetyp
309         do j=-maxthetyp,maxthetyp
310           do k=-maxthetyp,maxthetyp
311             aa0thet(i,j,k,iblock)=0.0d0
312             do l=1,ntheterm
313               aathet(l,i,j,k,iblock)=0.0d0
314             enddo
315             do l=1,ntheterm2
316               do m=1,nsingle
317                 bbthet(m,l,i,j,k,iblock)=0.0d0
318                 ccthet(m,l,i,j,k,iblock)=0.0d0
319                 ddthet(m,l,i,j,k,iblock)=0.0d0
320                 eethet(m,l,i,j,k,iblock)=0.0d0
321               enddo
322             enddo
323             do l=1,ntheterm3
324               do m=1,ndouble
325                 do mm=1,ndouble
326                  ffthet(mm,m,l,i,j,k,iblock)=0.0d0
327                  ggthet(mm,m,l,i,j,k,iblock)=0.0d0
328                 enddo
329               enddo
330             enddo
331           enddo
332         enddo
333       enddo
334       enddo
335       do iblock=1,2
336       do i=0,nthetyp
337         do j=-nthetyp,nthetyp
338           do k=-nthetyp,nthetyp
339             read (ithep,'(6a)') res1
340             read (ithep,*) aa0thet(i,j,k,iblock)
341             read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
342             read (ithep,*)
343      &       ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
344      &        (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
345      &        (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
346      &        (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle)
347      &        ,ll=1,ntheterm2)
348             read (ithep,*)
349      &      (((ffthet(llll,lll,ll,i,j,k,iblock),
350      &      ffthet(lll,llll,ll,i,j,k,iblock),
351      &         ggthet(llll,lll,ll,i,j,k,iblock)
352      &        ,ggthet(lll,llll,ll,i,j,k,iblock),
353      &         llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
354           enddo
355         enddo
356       enddo
357 C
358 C For dummy ends assign glycine-type coefficients of theta-only terms; the
359 C coefficients of theta-and-gamma-dependent terms are zero.
360 C
361       do i=1,nthetyp
362         do j=1,nthetyp
363           do l=1,ntheterm
364             aathet(l,i,j,nthetyp+1,iblock)=0.0d0
365             aathet(l,nthetyp+1,i,j,iblock)=0.0d0
366           enddo
367           aa0thet(i,j,nthetyp+1,iblock)=0.0d0
368           aa0thet(nthetyp+1,i,j,iblock)=0.0d0
369         enddo
370         do l=1,ntheterm
371           aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
372         enddo
373         aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
374       enddo
375       enddo
376 C Substitution for D aminoacids from symmetry.
377       do iblock=1,2
378       do i=-nthetyp,0
379         do j=-nthetyp,nthetyp
380           do k=-nthetyp,nthetyp
381            aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
382            do l=1,ntheterm
383            aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
384            enddo
385            do ll=1,ntheterm2
386             do lll=1,nsingle
387             bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
388             ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
389             ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
390             eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
391             enddo
392           enddo
393           do ll=1,ntheterm3
394            do lll=2,ndouble
395             do llll=1,lll-1
396             ffthet(llll,lll,ll,i,j,k,iblock)=
397      &      ffthet(llll,lll,ll,-i,-j,-k,iblock)
398             ffthet(lll,llll,ll,i,j,k,iblock)=
399      &      ffthet(lll,llll,ll,-i,-j,-k,iblock)
400             ggthet(llll,lll,ll,i,j,k,iblock)=
401      &      -ggthet(llll,lll,ll,-i,-j,-k,iblock)
402             ggthet(lll,llll,ll,i,j,k,iblock)=
403      &      -ggthet(lll,llll,ll,-i,-j,-k,iblock)
404             enddo !ll
405            enddo  !lll  
406           enddo   !llll
407          enddo    !k
408         enddo     !j
409        enddo      !i
410       enddo       !iblock
411
412 C
413 C Control printout of the coefficients of virtual-bond-angle potentials
414 C
415       if (lprint) then
416         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
417         do i=1,nthetyp+1
418           do j=1,nthetyp+1
419             do k=1,nthetyp+1
420               write (iout,'(//4a)')
421      &         'Type ',onelett(i),onelett(j),onelett(k)
422               write (iout,'(//a,10x,a)') " l","a[l]"
423               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
424               write (iout,'(i2,1pe15.5)')
425      &           (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
426             do l=1,ntheterm2
427               write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
428      &          "b",l,"c",l,"d",l,"e",l
429               do m=1,nsingle
430                 write (iout,'(i2,4(1pe15.5))') m,
431      &          bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
432      &          ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
433               enddo
434             enddo
435             do l=1,ntheterm3
436               write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
437      &          "f+",l,"f-",l,"g+",l,"g-",l
438               do m=2,ndouble
439                 do n=1,m-1
440                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
441      &              ffthet(n,m,l,i,j,k,iblock),
442      &              ffthet(m,n,l,i,j,k,iblock),
443      &              ggthet(n,m,l,i,j,k,iblock),
444      &              ggthet(m,n,l,i,j,k,iblock)
445                 enddo
446               enddo
447             enddo
448           enddo
449         enddo
450       enddo
451       call flush(iout)
452       endif
453 #endif
454
455 #ifdef CRYST_SC
456 C
457 C Read the parameters of the probability distribution/energy expression
458 C of the side chains.
459 C
460       do i=1,ntyp
461 cc      write (iout,*) "tu dochodze",i
462         read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
463         if (i.eq.10) then
464           dsc_inv(i)=0.0D0
465         else
466           dsc_inv(i)=1.0D0/dsc(i)
467         endif
468         if (i.ne.10) then
469         do j=1,nlob(i)
470           do k=1,3
471             do l=1,3
472               blower(l,k,j)=0.0D0
473             enddo
474           enddo
475         enddo  
476         bsc(1,i)=0.0D0
477         read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
478         censc(1,1,-i)=censc(1,1,i)
479         censc(2,1,-i)=censc(2,1,i)
480         censc(3,1,-i)=-censc(3,1,i)
481         do j=2,nlob(i)
482           read (irotam,*) bsc(j,i)
483           read (irotam,*) (censc(k,j,i),k=1,3),
484      &                                 ((blower(k,l,j),l=1,k),k=1,3)
485         censc(1,j,-i)=censc(1,j,i)
486         censc(2,j,-i)=censc(2,j,i)
487         censc(3,j,-i)=-censc(3,j,i)
488 C BSC is amplitude of Gaussian
489         enddo
490         do j=1,nlob(i)
491           do k=1,3
492             do l=1,k
493               akl=0.0D0
494               do m=1,3
495                 akl=akl+blower(k,m,j)*blower(l,m,j)
496               enddo
497               gaussc(k,l,j,i)=akl
498               gaussc(l,k,j,i)=akl
499              if (((k.eq.3).and.(l.ne.3))
500      &        .or.((l.eq.3).and.(k.ne.3))) then
501                 gaussc(k,l,j,-i)=-akl
502                 gaussc(l,k,j,-i)=-akl
503               else
504                 gaussc(k,l,j,-i)=akl
505                 gaussc(l,k,j,-i)=akl
506               endif
507             enddo
508           enddo 
509         enddo
510         endif
511       enddo
512       close (irotam)
513       if (lprint) then
514         write (iout,'(/a)') 'Parameters of side-chain local geometry'
515         do i=1,ntyp
516           nlobi=nlob(i)
517           if (nlobi.gt.0) then
518           write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
519      &     ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
520 c          write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
521 c          write (iout,'(a,f10.4,4(16x,f10.4))')
522 c     &                             'Center  ',(bsc(j,i),j=1,nlobi)
523 c          write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
524            write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
525      &                             'log h',(bsc(j,i),j=1,nlobi)
526            write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
527      &    'x',((censc(k,j,i),k=1,3),j=1,nlobi)
528 c          write (iout,'(a)')
529 c         do j=1,nlobi
530 c           ind=0
531 c           do k=1,3
532 c             do l=1,k
533 c              ind=ind+1
534 c              blower(k,l,j)=gaussc(ind,j,i)
535 c             enddo
536 c           enddo
537 c         enddo
538           do k=1,3
539             write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
540      &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
541           enddo
542           endif
543         enddo
544       endif
545 #else
546 C
547 C Read scrot parameters for potentials determined from all-atom AM1 calculations
548 C added by Urszula Kozlowska 07/11/2007
549 C
550       do i=1,ntyp
551         read (irotam,*)
552        if (i.eq.10) then
553          read (irotam,*)
554        else
555          do j=1,65
556            read(irotam,*) sc_parmin(j,i)
557          enddo
558        endif
559       enddo
560 #endif
561       close(irotam)
562 #ifdef CRYST_TOR
563 C
564 C Read torsional parameters in old format
565 C
566       read (itorp,*) ntortyp,nterm_old
567       write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
568       read (itorp,*) (itortyp(i),i=1,ntyp)
569       do i=1,ntortyp
570         do j=1,ntortyp
571           read (itorp,'(a)')
572           do k=1,nterm_old
573             read (itorp,*) kk,v1(k,j,i),v2(k,j,i) 
574           enddo
575         enddo
576       enddo
577       close (itorp)
578       if (lprint) then
579         write (iout,'(/a/)') 'Torsional constants:'
580         do i=1,ntortyp
581           do j=1,ntortyp
582             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
583             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
584           enddo
585         enddo
586       endif
587
588
589 #else
590 C
591 C Read torsional parameters
592 C
593       read (itorp,*) ntortyp
594       read (itorp,*) (itortyp(i),i=1,ntyp)
595       write (iout,*) 'ntortyp',ntortyp
596       do iblock=1,2
597       do i=-ntyp,-1
598        itortyp(i)=-itortyp(-i)
599       enddo
600 c      write (iout,*) 'ntortyp',ntortyp
601       do i=0,ntortyp-1
602         do j=-ntortyp+1,ntortyp-1
603           read (itorp,*) nterm(i,j,iblock),
604      &          nlor(i,j,iblock)
605           nterm(-i,-j,iblock)=nterm(i,j,iblock)
606           nlor(-i,-j,iblock)=nlor(i,j,iblock)
607           v0ij=0.0d0
608           si=-1.0d0
609           do k=1,nterm(i,j,iblock)
610             read (itorp,*) kk,v1(k,i,j,iblock),
611      &      v2(k,i,j,iblock)
612             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
613             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
614             v0ij=v0ij+si*v1(k,i,j,iblock)
615             si=-si
616          enddo
617           do k=1,nlor(i,j,iblock)
618             read (itorp,*) kk,vlor1(k,i,j),
619      &        vlor2(k,i,j),vlor3(k,i,j)
620             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
621           enddo
622           v0(i,j,iblock)=v0ij
623           v0(-i,-j,iblock)=v0ij
624         enddo
625       enddo
626       enddo
627       close (itorp)
628       if (lprint) then
629         write (iout,'(/a/)') 'Torsional constants:'
630         do i=1,ntortyp
631           do j=1,ntortyp
632             write (iout,*) 'ityp',i,' jtyp',j
633             write (iout,*) 'Fourier constants'
634             do k=1,nterm(i,j,iblock)
635               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
636      &        v2(k,i,j,iblock)
637             enddo
638             write (iout,*) 'Lorenz constants'
639             do k=1,nlor(i,j,iblock)
640               write (iout,'(3(1pe15.5))')
641      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
642             enddo
643           enddo
644         enddo
645       endif
646 C
647 C 6/23/01 Read parameters for double torsionals
648 C
649       do iblock=1,2
650       do i=0,ntortyp-1
651         do j=-ntortyp+1,ntortyp-1
652           do k=-ntortyp+1,ntortyp-1
653             read (itordp,'(3a1)') t1,t2,t3
654 c              write (iout,*) "OK onelett",
655 c     &         i,j,k,t1,t2,t3
656
657             if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
658      &        .or. t3.ne.toronelet(k)) then
659               write (iout,*) "Error in double torsional parameter file",
660      &         i,j,k,t1,t2,t3
661 #ifdef MPI
662               call MPI_Finalize(Ierror)
663 #endif
664                stop "Error in double torsional parameter file"
665             endif
666           read (itordp,*) ntermd_1(i,j,k,iblock),
667      &         ntermd_2(i,j,k,iblock)
668             ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
669             ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
670             read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,
671      &         ntermd_1(i,j,k,iblock))
672             read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,
673      &         ntermd_1(i,j,k,iblock))
674             read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,
675      &         ntermd_1(i,j,k,iblock))
676             read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,
677      &         ntermd_1(i,j,k,iblock))
678 C Martix of D parameters for one dimesional foureir series
679             do l=1,ntermd_1(i,j,k,iblock)
680              v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
681              v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
682              v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
683              v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
684 c            write(iout,*) "whcodze" ,
685 c     & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
686             enddo
687             read (itordp,*) ((v2c(l,m,i,j,k,iblock),
688      &         v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
689      &         v2s(m,l,i,j,k,iblock),
690      &         m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
691 C Martix of D parameters for two dimesional fourier series
692             do l=1,ntermd_2(i,j,k,iblock)
693              do m=1,l-1
694              v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
695              v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
696              v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
697              v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
698              enddo!m
699             enddo!l
700           enddo!k
701         enddo!j
702       enddo!i
703       enddo!iblock
704       if (lprint) then
705       write (iout,*)
706       write (iout,*) 'Constants for double torsionals'
707       do iblock=1,2
708       do i=0,ntortyp-1
709         do j=-ntortyp+1,ntortyp-1
710           do k=-ntortyp+1,ntortyp-1
711             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
712      &        ' nsingle',ntermd_1(i,j,k,iblock),
713      &        ' ndouble',ntermd_2(i,j,k,iblock)
714             write (iout,*)
715             write (iout,*) 'Single angles:'
716             do l=1,ntermd_1(i,j,k,iblock)
717               write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
718      &           v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
719      &           v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
720      &           v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
721             enddo
722             write (iout,*)
723             write (iout,*) 'Pairs of angles:'
724             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
725             do l=1,ntermd_2(i,j,k,iblock)
726               write (iout,'(i5,20f10.5)')
727      &         l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
728             enddo
729             write (iout,*)
730            write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
731             do l=1,ntermd_2(i,j,k,iblock)
732               write (iout,'(i5,20f10.5)')
733      &         l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
734      &         (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
735             enddo
736             write (iout,*)
737           enddo
738         enddo
739       enddo
740       enddo
741       endif
742 #endif
743 C Read of Side-chain backbone correlation parameters
744 C Modified 11 May 2012 by Adasko
745 CCC
746 C
747       read (isccor,*) nsccortyp
748       read (isccor,*) (isccortyp(i),i=1,ntyp)
749       do i=-ntyp,-1
750         isccortyp(i)=-isccortyp(-i)
751       enddo
752       iscprol=isccortyp(20)
753 c      write (iout,*) 'ntortyp',ntortyp
754       maxinter=3
755 cc maxinter is maximum interaction sites
756       do l=1,maxinter
757       do i=1,nsccortyp
758         do j=1,nsccortyp
759           read (isccor,*)
760      &nterm_sccor(i,j),nlor_sccor(i,j)
761           write (iout,*) nterm_sccor(i,j)
762           v0ijsccor=0.0d0
763           v0ijsccor1=0.0d0
764           v0ijsccor2=0.0d0
765           v0ijsccor3=0.0d0
766           si=-1.0d0
767           nterm_sccor(-i,j)=nterm_sccor(i,j)
768           nterm_sccor(-i,-j)=nterm_sccor(i,j)
769           nterm_sccor(i,-j)=nterm_sccor(i,j)
770           write (iout,*) nterm_sccor(i,j),nterm_sccor(-i,j),
771      &    nterm_sccor(-i,-j),nterm_sccor(i,-j)
772           do k=1,nterm_sccor(i,j)
773             read (isccor,*) kk,v1sccor(k,l,i,j)
774      &    ,v2sccor(k,l,i,j)
775             if (j.eq.iscprol) then
776              if (i.eq.isccortyp(10)) then
777              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
778              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
779              else
780              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
781      &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
782              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
783      &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
784              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
785              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
786              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
787              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
788              endif
789             else
790              if (i.eq.isccortyp(10)) then
791              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
792              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
793              else
794                if (j.eq.isccortyp(10)) then
795              v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
796              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
797                else
798              v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
799              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
800              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
801              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
802              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
803              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
804                 endif
805                endif
806             endif
807             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
808             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
809             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
810             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
811             si=-si
812            enddo
813           do k=1,nlor_sccor(i,j)
814             read (isccor,*) kk,vlor1sccor(k,i,j),
815      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
816             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
817      &(1+vlor3sccor(k,i,j)**2)
818           enddo
819           v0sccor(l,i,j)=v0ijsccor
820           v0sccor(l,-i,j)=v0ijsccor1
821           v0sccor(l,i,-j)=v0ijsccor2
822           v0sccor(l,-i,-j)=v0ijsccor3
823           enddo
824         enddo
825       enddo
826       close (isccor)
827       if (lprint) then
828         write (iout,'(/a/)') 'Torsional constants of SCCORR:'
829         do i=1,nsccortyp
830           do j=1,nsccortyp
831             write (iout,*) 'ityp',i,' jtyp',j
832             write (iout,*) 'Fourier constants'
833             do k=1,nterm_sccor(i,j)
834               write (iout,'(2(1pe15.5))')
835      & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
836             enddo
837             write (iout,*) 'Lorenz constants'
838             do k=1,nlor_sccor(i,j)
839               write (iout,'(3(1pe15.5))')
840      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
841             enddo
842           enddo
843         enddo
844       endif
845 C
846 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
847 C         interaction energy of the Gly, Ala, and Pro prototypes.
848 C
849       read (ifourier,*) nloctyp
850       do i=0,nloctyp-1
851         read (ifourier,*)
852         read (ifourier,*) (b(ii,i),ii=1,13)
853         if (lprint) then
854         write (iout,*) 'Type',i
855         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
856         endif
857         B1(1,i)  = b(3,i)
858         B1(2,i)  = b(5,i)
859         B1(1,-i) = b(3,i)
860         B1(2,-i) = -b(5,i)
861 c        b1(1,i)=0.0d0
862 c        b1(2,i)=0.0d0
863         B1tilde(1,i) = b(3,i)
864         B1tilde(2,i) =-b(5,i)
865         B1tilde(1,-i) =-b(3,i)
866         B1tilde(2,-i) =b(5,i)
867 c        b1tilde(1,i)=0.0d0
868 c        b1tilde(2,i)=0.0d0
869         B2(1,i)  = b(2,i)
870         B2(2,i)  = b(4,i)
871         B2(1,-i)  =b(2,i)
872         B2(2,-i)  =-b(4,i)
873
874 c        b2(1,i)=0.0d0
875 c        b2(2,i)=0.0d0
876         CC(1,1,i)= b(7,i)
877         CC(2,2,i)=-b(7,i)
878         CC(2,1,i)= b(9,i)
879         CC(1,2,i)= b(9,i)
880         CC(1,1,-i)= b(7,i)
881         CC(2,2,-i)=-b(7,i)
882         CC(2,1,-i)=-b(9,i)
883         CC(1,2,-i)=-b(9,i)
884 c        CC(1,1,i)=0.0d0
885 c        CC(2,2,i)=0.0d0
886 c        CC(2,1,i)=0.0d0
887 c        CC(1,2,i)=0.0d0
888         Ctilde(1,1,i)=b(7,i)
889         Ctilde(1,2,i)=b(9,i)
890         Ctilde(2,1,i)=-b(9,i)
891         Ctilde(2,2,i)=b(7,i)
892         Ctilde(1,1,-i)=b(7,i)
893         Ctilde(1,2,-i)=-b(9,i)
894         Ctilde(2,1,-i)=b(9,i)
895         Ctilde(2,2,-i)=b(7,i)
896
897 c        Ctilde(1,1,i)=0.0d0
898 c        Ctilde(1,2,i)=0.0d0
899 c        Ctilde(2,1,i)=0.0d0
900 c        Ctilde(2,2,i)=0.0d0
901         DD(1,1,i)= b(6,i)
902         DD(2,2,i)=-b(6,i)
903         DD(2,1,i)= b(8,i)
904         DD(1,2,i)= b(8,i)
905         DD(1,1,-i)= b(6,i)
906         DD(2,2,-i)=-b(6,i)
907         DD(2,1,-i)=-b(8,i)
908         DD(1,2,-i)=-b(8,i)
909 c        DD(1,1,i)=0.0d0
910 c        DD(2,2,i)=0.0d0
911 c        DD(2,1,i)=0.0d0
912 c        DD(1,2,i)=0.0d0
913         Dtilde(1,1,i)=b(6,i)
914         Dtilde(1,2,i)=b(8,i)
915         Dtilde(2,1,i)=-b(8,i)
916         Dtilde(2,2,i)=b(6,i)
917         Dtilde(1,1,-i)=b(6,i)
918         Dtilde(1,2,-i)=-b(8,i)
919         Dtilde(2,1,-i)=b(8,i)
920         Dtilde(2,2,-i)=b(6,i)
921
922 c        Dtilde(1,1,i)=0.0d0
923 c        Dtilde(1,2,i)=0.0d0
924 c        Dtilde(2,1,i)=0.0d0
925 c        Dtilde(2,2,i)=0.0d0
926         EE(1,1,i)= b(10,i)+b(11,i)
927         EE(2,2,i)=-b(10,i)+b(11,i)
928         EE(2,1,i)= b(12,i)-b(13,i)
929         EE(1,2,i)= b(12,i)+b(13,i)
930         EE(1,1,-i)= b(10,i)+b(11,i)
931         EE(2,2,-i)=-b(10,i)+b(11,i)
932         EE(2,1,-i)=-b(12,i)+b(13,i)
933         EE(1,2,-i)=-b(12,i)-b(13,i)
934
935 c        ee(1,1,i)=1.0d0
936 c        ee(2,2,i)=1.0d0
937 c        ee(2,1,i)=0.0d0
938 c        ee(1,2,i)=0.0d0
939 c        ee(2,1,i)=ee(1,2,i)
940
941       enddo
942       if (lprint) then
943       do i=1,nloctyp
944         write (iout,*) 'Type',i
945         write (iout,*) 'B1'
946 c        write (iout,'(f10.5)') B1(:,i)
947         write(iout,*) B1(1,i),B1(2,i)
948         write (iout,*) 'B2'
949 c        write (iout,'(f10.5)') B2(:,i)
950         write(iout,*) B2(1,i),B2(2,i)
951         write (iout,*) 'CC'
952         do j=1,2
953           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
954         enddo
955         write(iout,*) 'DD'
956         do j=1,2
957           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
958         enddo
959         write(iout,*) 'EE'
960         do j=1,2
961           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
962         enddo
963       enddo
964       endif
965
966 C Read electrostatic-interaction parameters
967 C
968       if (lprint) then
969         write (iout,'(/a)') 'Electrostatic interaction constants:'
970         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') 
971      &            'IT','JT','APP','BPP','AEL6','AEL3'
972       endif
973       read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
974       read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
975       read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
976       read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
977       close (ielep)
978       do i=1,2
979         do j=1,2
980         rri=rpp(i,j)**6
981         app (i,j)=epp(i,j)*rri*rri 
982         bpp (i,j)=-2.0D0*epp(i,j)*rri
983         ael6(i,j)=elpp6(i,j)*4.2D0**6
984         ael3(i,j)=elpp3(i,j)*4.2D0**3
985         lprint=.true.
986         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
987      &                    ael6(i,j),ael3(i,j)
988          lprint=.false.
989         enddo
990       enddo
991 C
992 C Read side-chain interaction parameters.
993 C
994       read (isidep,*) ipot,expon
995       if (ipot.lt.1 .or. ipot.gt.5) then
996         write (iout,'(2a)') 'Error while reading SC interaction',
997      &               'potential file - unknown potential type.'
998         stop
999       endif
1000       expon2=expon/2
1001       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1002      & ', exponents are ',expon,2*expon 
1003       goto (10,20,30,30,40) ipot
1004 C----------------------- LJ potential ---------------------------------
1005    10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1006       if (lprint) then
1007         write (iout,'(/a/)') 'Parameters of the LJ potential:'
1008         write (iout,'(a/)') 'The epsilon array:'
1009         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1010         write (iout,'(/a)') 'One-body parameters:'
1011         write (iout,'(a,4x,a)') 'residue','sigma'
1012         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1013       endif
1014       goto 50
1015 C----------------------- LJK potential --------------------------------
1016    20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1017      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1018       if (lprint) then
1019         write (iout,'(/a/)') 'Parameters of the LJK potential:'
1020         write (iout,'(a/)') 'The epsilon array:'
1021         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1022         write (iout,'(/a)') 'One-body parameters:'
1023         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
1024         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1025      &        i=1,ntyp)
1026       endif
1027       goto 50
1028 C---------------------- GB or BP potential -----------------------------
1029    30 do i=1,ntyp
1030        read (isidep,*)(eps(i,j),j=i,ntyp)
1031       enddo
1032       read (isidep,*)(sigma0(i),i=1,ntyp)
1033       read (isidep,*)(sigii(i),i=1,ntyp)
1034       read (isidep,*)(chip(i),i=1,ntyp)
1035       read (isidep,*)(alp(i),i=1,ntyp)
1036 C For the GB potential convert sigma'**2 into chi'
1037       if (ipot.eq.4) then
1038         do i=1,ntyp
1039           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1040         enddo
1041       endif
1042       if (lprint) then
1043         write (iout,'(/a/)') 'Parameters of the BP potential:'
1044         write (iout,'(a/)') 'The epsilon array:'
1045         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1046         write (iout,'(/a)') 'One-body parameters:'
1047         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
1048      &       '    chip  ','    alph  '
1049         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1050      &                     chip(i),alp(i),i=1,ntyp)
1051       endif
1052       goto 50
1053 C--------------------- GBV potential -----------------------------------
1054    40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1055      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1056      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1057       if (lprint) then
1058         write (iout,'(/a/)') 'Parameters of the GBV potential:'
1059         write (iout,'(a/)') 'The epsilon array:'
1060         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1061         write (iout,'(/a)') 'One-body parameters:'
1062         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
1063      &      's||/s_|_^2','    chip  ','    alph  '
1064         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1065      &           sigii(i),chip(i),alp(i),i=1,ntyp)
1066       endif
1067    50 continue
1068       close (isidep)
1069 C-----------------------------------------------------------------------
1070 C Calculate the "working" parameters of SC interactions.
1071       do i=2,ntyp
1072         do j=1,i-1
1073           eps(i,j)=eps(j,i)
1074         enddo
1075       enddo
1076       do i=1,ntyp
1077         do j=i,ntyp
1078           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1079           sigma(j,i)=sigma(i,j)
1080           rs0(i,j)=dwa16*sigma(i,j)
1081           rs0(j,i)=rs0(i,j)
1082         enddo
1083       enddo
1084       if (lprint) write (iout,'(/a/10x,7a/72(1h-))') 
1085      & 'Working parameters of the SC interactions:',
1086      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
1087      & '  chi1   ','   chi2   ' 
1088       do i=1,ntyp
1089         do j=i,ntyp
1090           epsij=eps(i,j)
1091           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1092             rrij=sigma(i,j)
1093           else
1094             rrij=rr0(i)+rr0(j)
1095           endif
1096           r0(i,j)=rrij
1097           r0(j,i)=rrij
1098           rrij=rrij**expon
1099           epsij=eps(i,j)
1100           sigeps=dsign(1.0D0,epsij)
1101           epsij=dabs(epsij)
1102           aa(i,j)=epsij*rrij*rrij
1103           bb(i,j)=-sigeps*epsij*rrij
1104           aa(j,i)=aa(i,j)
1105           bb(j,i)=bb(i,j)
1106           if (ipot.gt.2) then
1107             sigt1sq=sigma0(i)**2
1108             sigt2sq=sigma0(j)**2
1109             sigii1=sigii(i)
1110             sigii2=sigii(j)
1111             ratsig1=sigt2sq/sigt1sq
1112             ratsig2=1.0D0/ratsig1
1113             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1114             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1115             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1116           else
1117             rsum_max=sigma(i,j)
1118           endif
1119 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1120             sigmaii(i,j)=rsum_max
1121             sigmaii(j,i)=rsum_max 
1122 c         else
1123 c           sigmaii(i,j)=r0(i,j)
1124 c           sigmaii(j,i)=r0(i,j)
1125 c         endif
1126 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1127           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1128             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1129             augm(i,j)=epsij*r_augm**(2*expon)
1130 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1131             augm(j,i)=augm(i,j)
1132           else
1133             augm(i,j)=0.0D0
1134             augm(j,i)=0.0D0
1135           endif
1136           if (lprint) then
1137             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') 
1138      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1139      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1140           endif
1141         enddo
1142       enddo
1143 C
1144 C Define the SC-p interaction constants
1145 C
1146 #ifdef OLDSCP
1147       do i=1,20
1148 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates 
1149 C helix formation)
1150 c       aad(i,1)=0.3D0*4.0D0**12
1151 C Following line for constants currently implemented
1152 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1153         aad(i,1)=1.5D0*4.0D0**12
1154 c       aad(i,1)=0.17D0*5.6D0**12
1155         aad(i,2)=aad(i,1)
1156 C "Soft" SC-p repulsion
1157         bad(i,1)=0.0D0
1158 C Following line for constants currently implemented
1159 c       aad(i,1)=0.3D0*4.0D0**6
1160 C "Hard" SC-p repulsion
1161         bad(i,1)=3.0D0*4.0D0**6
1162 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1163         bad(i,2)=bad(i,1)
1164 c       aad(i,1)=0.0D0
1165 c       aad(i,2)=0.0D0
1166 c       bad(i,1)=1228.8D0
1167 c       bad(i,2)=1228.8D0
1168       enddo
1169 #else
1170 C
1171 C 8/9/01 Read the SC-p interaction constants from file
1172 C
1173       do i=1,ntyp
1174         read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1175       enddo
1176       do i=1,ntyp
1177         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1178         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1179         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1180         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1181       enddo
1182
1183       if (lprint) then
1184         write (iout,*) "Parameters of SC-p interactions:"
1185         do i=1,20
1186           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1187      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1188         enddo
1189       endif
1190 #endif
1191 C
1192 C Define the constants of the disulfide bridge
1193 C
1194       ebr=-5.50D0
1195 c
1196 c Old arbitrary potential - commented out.
1197 c
1198 c      dbr= 4.20D0
1199 c      fbr= 3.30D0
1200 c
1201 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1202 c energy surface of diethyl disulfide.
1203 c A. Liwo and U. Kozlowska, 11/24/03
1204 c
1205       D0CM = 3.78d0
1206       AKCM = 15.1d0
1207       AKTH = 11.0d0
1208       AKCT = 12.0d0
1209       V1SS =-1.08d0
1210       V2SS = 7.61d0
1211       V3SS = 13.7d0
1212
1213       if (lprint) then
1214       write (iout,'(/a)') "Disulfide bridge parameters:"
1215       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1216       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1217       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1218       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1219      & ' v3ss:',v3ss
1220       endif
1221       return
1222       end