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